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Introduction 

High-intensity  focused  ultrasound  (H1FU)  is  gaining  wider  acceptance  in  noninvasive  or  minimally 
invasive  targeting  of  abnormal  tissues  (e.g.  cancer)  for  destruction.  Piezocomposite  transducer 
technology,  especially  for  phased  arrays,  is  providing  high-quality  HIFU  applicators  with  increased 
bandwidth  and  reduced  parasitic  cross  coupling  between  the  array  elements.  In  addition  to  increasing 
the  efficacy  of  HIFU  applicators,  these  technological  enhancements  allow  for  the  use  of  HIFU  arrays 
in  imaging  the  target  region  before,  after,  and  intermittently  during  lesion  formation.  This  leads  to  a 
unique  paradigm  of  image-guided  surgery  with  HIFU  in  which  the  coordinate  systems  for  both 
therapy  and  imaging  are  inherently  registered.  This  project  investigates  the  feasibility  of  using 
piezocomposite  phased  arrays  as  dual-mode  applicators  for  the  noninvasive  treatment  of  primary 
breast  cancers.  Both  therapeutic  and  imaging  capabilities  of  the  dual-mode  arrays  are  investigated 
leading  to  a  real-time  dual-mode  array  system  to  be  used  in  pursuing  in  vivo  animal  experiments  in 
the  future. 

Body 

This  report  is  structured  in  accordance  with  the  approved  statement  of  work  (SoW).  In  what  follows, 
we  give  the  tasks  and  subtasks  of  the  approved  SoW  with  each  subtask  followed  directly  by  what  has 
been  accomplished  with  respect  to  it  in  year  2.  For  the  subtasks  planned  for  years  1  and  3  of  the 
grant  period,  they  are  given  here  for  completeness,  with  short  comments  indicating  their  status,  c.g. 
sublasks  completed  in  year  1  as  indicated  in  year  1  progress  report. 

Task  1.  Thresholds  for  Thermal  Ablation  of  Breast  and  Fatty  Tissue  (Months  1  12): 

a)  Investigation  of  the  intensity/exposure  threshold  curve;  lesion  size  and  characterization 
of  damage  (1  -  3):  Completed  (Year  1  Report) 

b)  Imaging  of  discrete  thermal  lesions  with  therapeutic  arrays  (1  -  3):  Completed  (Year  1 
Report). 

c)  Long-duration  volumetric  ablation  of  porcine  fatty  tissue  (3  -  9):  Completed  (Year  1 
Report). 

d)  Imaging  of  volumetric  ablations  using  the  therapeutic  array  and  diagnostic  scanners  (3 
-  12):  Completed  (Appendix  IV). 

Task  2.  Treatment  Planning  and  Optimization  of  Volumetric  Ablation  with  Phased  Arrays  in  Fatly 
Tissue  (Months  6  -  24): 

a)  Thermal  modeling  based  on  bioheat  equation  and  Saparato  and  Dewey ■  thermal  dose 
integral  for  damage:  Discrete  lesions  (6-12):  Completed  (Year  1  Report) 

b)  Thermal  modeling  for  multiple  lesions  with  variable  levels  of  proximity  and  cooling 
time  between  shots  (12  -  18):  We  have  simulated  several  treatment  scenarios  feasible 
with  the  dual  mode  arrays.  Starting  with  the  conventional  procedure  of  forming  a 
volumetric  lesion  by  doing  a  series  of  discrete  lesions  one  at  a  time  with  sufficient  wait 
time  between  lesions  to  allow  for  the  temperature  to  cool  down  (close  to  base  line). 
This  technique  requires  long  waiting  time  between  (thus  leading  to  elongating  the 
overall  treatment  time)  if  one  is  to  avoid  excessive  heating  in  the  prefocal  region, 
which  could  lead  to  unintended  damage  of  intervening  tissue.  We  have  also 
investigated  the  newly  emerging  technique  of  (electronically)  raster  scanning  the 
therapeutic  beam  continuously  at  speeds  of  1  -  3  mm/s  to  form  volumetric  lesion. 
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While  the  origin  of  this  technique  was  proposed  by  Dr.  ter  Haar  and  first  used 
clinically  on  the  Chongqing  HIFU  system,  it  was  very  nicely  illustrated  in  vitro  by  the 
group  at  the  University  of  Washington  in  their  recent  presentations.  We  have  also 
carried  out  some  tissue  heating  experiments  performing  the  raster  scan  approach  to 
confirm  the  simulation  results.  The  results  confirm  that  the  raster  scan  approach  is 

superior  to  the  shoot-and-wait  approach  in  terms  of  treatment  time  and  prcfocal 
heating. 

c)  Optimization  of  multiple-focus  phased  array  patterns  for  simultaneous  placement 
multiple  discrete  lesions  (12  -  24):  We  have  also  begun  investigating  the  use  of 
multiple  focus  patterns  for  the  minimization  of  collateral  damage  in  the  prcfocal 
region.  This  follows  work  done  much  earlier  by  the  PI  in  investigating  phased  arrays 
for  hyperthermia  [Ebbini:91].  For  the  64-element  prototype  used  in  this  investigation, 
multiple-focus  patterns  do  not  give  a  distinct  advantage  over  the  raster  scan  method 
described  in  the  previous  subtask.  This  is  not  surprising  given  the  relatively  small 
number  of  array  elements  allowing  for  only  a  small  number  of  simultaneous  foci  (2 
4).  Simulations  of  other  2D  arrays  with  larger  number  of  elements,  however,  continue 
to  show  that  multiple-focus  patterns  offer  significant  advantage  in  reducing  prefocal 
heating  compared  with  any  form  of  mechanical  scanning  (including  the  raster  scan 
method). 

Task  3.  Detection  and  Localization  of  Cavitation  Activity  During  Thermal  Lesion  Formation  in 
Breast  and  Fatty  Tissue  (Months  1-24): 

a)  Detection  of  subharmonic  activity  in  single-channel  and  heamformed  data  (I  -  12): 
Completed  (Year  1  Report). 

b)  Localization  of  cavitation  from  beamforming  of  multiple  receiving  channels  (12  —  24): 
We  are  still  unable  to  demonstrate  the  sensitivity  of  the  array  elements  to  subharmonic 
generation  using  short  (microsecond)  imaging  pulses  suitable  for  localization. 
However,  we  still  think  that  localization  of  cavitation  events  using  microsecond  pulses 
is  feasible.  We  intend  to  order  custom  designed  transducers,  similar  to  the  prototype  we 
already  have,  with  sufficient  number  of  hydrophones  integrated  with  the  transducer. 
The  new  design  will  attempt  to  improve  the  isolation  between  the  therapeutic  array 
elements  and  the  receiving  hydrophones  (so  they  can  be  used  simultaneously  with  the 
therapeutic  excitation).  Furthermore,  the  hydrophones  will  be  designed  to  have 
maximum  sensitivity  at  the  subharmonic.  The  use  of  an  array  of  receiving  hydrophones 
(on  the  order  of  8  elements)  will  still  allow  us  to  achieve  localization  (provided  that  wc 
can  consistently  produce  cavitational  oscillations  with  microsecond  excitation  pulses). 

c)  Localization  using  time-frequency  and  related  methods  (6  -  18):  Our  intuition  with 

regard  to  the  harmonic  content  of  echoes  from  HIFU-induced  thermal  lesions  is 
supported  by  data  from  numerous  lesion  formation  experiments.  We  have  believed  all 
along  that  HIFU  beams  produce  gas  bubbles  in  the  target  tissue  within  the  beam  waist 
in  the  focal  region.  At  normal  exposure,  these  are  stable  microbubbles  that  can  be 
beneficial  in  enhancing  the  power  deposition  at  the  target.  They  are  also  responsible  for 
significant  increase  in  echogenecity  and  second  harmonic  generation,  especially  when 
imaged  at  the  therapeutic  frequency.  The  figure  below  shows  25  dB  grayscale  images 
of  ex  vivo  liver  tissue  before  and  after  lesion  formation  using  a  3  second  exposure  at 
850  W/cnr.  One  can  see  the  increase  in  echogenecity  at  the  lesion  location  (92  102 
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mm  axial).  On  the  right  hand  side,  spectrograms  from  the  RF  line  passing  through  the 
lesion  (indicated  by  the  arrows)  are  computed  with  60  dB  dynamic  range.  The 
spectrogram  produced  from  the  image  after  lesion  formation  shows  high  level  of 
second  harmonic  activity  localized  to  the  lesion  boundaries.  It  is  this  kind  of 
observation  that  led  to  the  development  of  the  second  harmonic  imaging  and. 
subsequently,  to  the  quadratic  imaging  using  the  second  order  Volterra  filters. 


Task  4.  Image  Characterization  of  Thermal  Lesions  in  Breast  and  Fatty  Tissue  (Months  1  -  24) 

a)  Characterization  of  grayscale  images  for  discrete  thermal  lesions  (l  -  6):  Completed  (Year  1 
Report). 

b)  Characterization  of  tissue  dependent  parameters  (I  -  12):  Completed  (Year  1  Report: 
(Expected)  negative  results). 

c)  Characterization  of  second  harmonic  imaging  (12  -  24):  Over  the  past  18  months,  we  have 

focused  our  attention  on  quadratic  imaging  which  results  from  the  second  order  Volterra  filter 
formulation.  The  application  of  this  filter  to  the  lesion  visualization  is  described  in  Appendix  II 
and  Appendix  IV  for  image  data  acquired  with  a  commercial  diagnostic  system.  Appendix  III 
describe  its  application  to  image  data  acquired  with  the  dual-mode  array.  A  fuller  description  is 
provided  in  the  patent  application  and  the  journal  paper.  The 

basic  idea  behind  quadratic  filtering  is  illustrated  by  the  figure  below  whereby  the  beamformed 
Rf-  data  from  the  array,  y(n),  is  decomposed  into  a  linear  and  quadratic  components  as  shown. 
We  have  developed  an  algorithm  for  deriving  the  linear  and  quadratic  filters,  hi,  and  h0 .  shown 
in  the  figure. 
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d)  C™'iuWn  °f  imaflS  Parameters  with  histologic  characterization  of  tissue  ablations  (12  - 
"  '  hlS  ,s  current,y  on§olng  work.  We  have  requested  the  addition  of  Dr.  Robert  Griffin  from 
the  Department  of  Therapeutic  Radiology  as  an  investigator  on  this  project  to  help  with  this 

task  beginning  the  3  year  of  the  project.  We  have  not  yet  heard  from  the  program  manager 
with  regards  to  this  request. 


(Morahs  !TS24b)aSed  AdaP'iVe  Ref0CUSlng  0f  TheraPculic  Arrays  ™  Inhomogeneous  Fatty  Tissue 

Focusing  with  reference  to  natural  specular  targets  (12  -  18): 

1  he  idea  behind  this  task  may  be  illustrated  easily  by  Figure  1,  which 
shows  grayscale  (25  dB)  images  of  a  fresh  (untreated)  sample  of 
porcine  liver.  The  image  on  the  left  hand  side  was  acquired  using 
smgle-transmit  focus  at  the  geometric  center.  Two  large  reflectors 
(blood  vessels)  and  two  smaller  ones  can  be  seen  in  the  image.  Using 
this  image,  the  coordinates  of  the  reflector  at  (5,90)  mm  were  obtained. 

I  his  was  done  through  an  interface  program  running  under  MATLAB 
v\  hereby  the  user  clicks  the  mouse  on  a  target  point  in  the  image  and 
the  program  finds  the  coordinates  of  the  target  based  on  the 
beamiorming  data  that  produced  the  image  in  the  first  place.  We  have 
then  used  the  coordinates  of  the  target  to  compute  the  delays  needed  to 
locus  the  transmit  beam  on  this  target.  These  delays  were  then  used  to 
produce  the  single  transmit  beam  needed  to  obtain  the  image  on  the 
right  hand  side  (25  dB  dynamic  range).  The  image  clearly  shows  that 
the  target  is  illuminated  while  the  other  three  vessels  are  relatively 

dim.  especially  the  vessel  at  (-5,90+).  This  result  is  quite  significant  as 
ii  show's  that: 

1 .  Critical  structures  (e.g.,  significant  blood  vessels)  in  the  vicinity  of  the  focus  point  of  the  transmit  focus  are 
i  uminated  by  its  sidelobe  and  can  be  seen  in  the  single-transmit  imaging  mode. 

-  Spatial  coordinates  of  these  critical  structures  can  be  accurately  estimated  from  the  beamformed  single-transmil 


Figure  1  The  image  to  the  left  is 
obtained  using  a  single  transmit  focus 
along  the  central  axis.  The  image  to  the 
right  is  obtained  by  a  single  transmit 
focus  at  the  blood  vessel  at  90  mm  axial 
and  5  mm  lateral. 


3.  Based  on  their  estimated  coordinates,  these  critical  structures  can  be  targeted  (or  avoided)  by  synthesizin'- 

transmit  beams  with  maxima  (or  nulls)  at  these  points.  & 

Figure  1  shows  the  case  of  targeting  a  blood  vessel  with  the  transmit  beam.  One  could  have  just  as  easily  generated 
a  earn  that  generates  a  minimum  at  the  vessel  location  (avoidance).  Many  other  scenarios  can  be  envisioned  once 

this  capability  is  developed  for  guidance. 

Temperature  feedback  (12  -  24):  The  sensitivity  of  the  speckle  component  of  the  dual-mode  array  to 
temperature  changes  has  been  established  experimentally.  We  have  developed  two  different 
temperature  imaging  algorithms  (based  on  pulse-echo  ultrasound)  under  previous  NIH  funding  (Seip 

and  Ebbini  and  Simon  and  Ebbini).  We  fully  expect  that  temperature  feedback  will  be  implemented 
with  the  dual-mode  array  system. 


iscrete  thermal  lesions  as  beacons  (12  -  24):  This  works  in  a  similar  manner  to  the  illustration  in 
figure  1.  However,  it  will  not  work  with  extended  lesions.  Our  original  plan  was  to  generate  a  very 
small  lesion  that  appears  as  a  point  reflector  in  the  image  and  use  it  as  a  reference  point.  Our  results  so 
iar  provide  a  proof  of  concept.  It  is  expected  that  the  actual  implementation  of  this  procedure  will  be 
tied  closely  to  establishing  the  thresholds  for  damage  in  the  target  tissue,  e.g.  breast  tissue. 


Task  6-  Real-time  Dual-mode  Phased  Array  System  for  Volumetric  Thermal  Ablation  of  Breast  and 
Fatty  Tissue  (Months  1  -  36) 
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a)  Design  and  fabrication  of  64-channel  receive  system  (1  -  12):  A  new  receiver  circuit 
employing  multiplexers  for  minimizing  the  number  of  analog-to-digital  (A/D)  converters  has 
been  designed  and  fabricated.  This  will  significantly  reduce  the  cost  and  complexity  of  the 
receive  electronics.  High-speed  multiplexers  allowed  us  to  perform  8:1  or  4:1  multiplexing  to 
each  A/D  while  maintaining  a  minimum  sampling  rate  of  10  MHz.  This  is  sufficient  for  the 
dual-mode  array  operating  at  1  MHz  (and  can  capture  the  second  harmonic  without  aliasing). 
Results  from  this  new  circuit  will  be  reported  at  the  Ultrasonics  Symposium  in  October. 

b)  Design  and  fabrication  of  transmit/receive  control  circuitry  (1  18):  Completed  (Year  1 

Report). 

c)  DSP-based  real-time  beamformer  (12  -  30):  We  have  identified  a  VXI  based  real-time 
controller  for  beamforming.  This  controller  will  be  sufficiently  fast  to  perform  real-time 
beamforming  and  uploading  of  image  data  to  the  controller  workstation. 

d)  Experimental  testing,  calibration,  and  characterization  of  imaging  system  (24  -  36):  We  are 
on  track  for  testing  and  characterization  of  the  system  for  real-time  data  acquisition  from  all 
channels  during  the  first  half  of  the  third  year.  We  expect  to  be  able  to  form  images  with  a 
frame  rate  on  the  order  of  1 0  -  1 5  frame  per  second.  If  this  is  achieved,  we  will  be  in  an 
excellent  shape  to  propose  the  use  of  the  dual-mode  system  in  vivo  animal  experiments  as 
part  of  a  future  grant  (possibly  from  the  NIH). 


Key  Research  Accomplishments 

•  We  have  established  that  the  real-time  harmonic  activity  during  and  after  lesion  formation  is 
observable  using  diagnostic  scanners  (Technos  MPX,  Esaote,  S.p.A.,  Genoa  Italy).  Hypothesis 
on  microbubble  activity  is  strengthened  by  real-time  imaging  results  using  the  pulse  inversion 
technique  (Appendix  II). 

•  We  have  extended  the  second  harmonic  imaging  results  to  quadratic  imaging  using  the  second- 
order  Volterra  filter.  This  nonlinear  imaging  mode  is  superior  to  second  harmonic  imaging  in 
that  it  improves  spatial  and  contrast  resolution  and  virtually  eliminates  beamforming  artifacts. 
The  latter  severely  limits  the  dynamic  range  of  the  2nd  harmonic  imaging  technique  in  single- 
transmit  imaging  mode.  (Appendix  II,  III,  IV). 

•  We  have  established  the  correlation  between  pulse  inversion  data  and  quadratic  data  for  discrete 
lesions  (Appendix  II)  and  volumetric  lesions  (Appendix  IV). 

•  We  have  demonstrated  that  nonlinear  imaging  methods  like  pulse  inversion  and  quadratic 
filtering  produce  more  accurate  maps  of  the  actual  lesions  (Appendix  IV). 

•  Full  characterization  of  the  imaging  capabilities  of  our  64-element  prototype  dual-mode  array. 
Comparison  with  diagnostic  imaging  using  the  Esaote  Technos  MPX  scanner  (CA421  convex  for 
abdominal  imaging).  This  will  be  reported  at  the  2003  International  Ultrasonics  Symposium  in 
October. 

•  Characterization  of  the  transient  behavior  of  the  harmonic  contents  of  the  echo  signal  from  the 
lesion  location  in  ex  vivo  tissue  (several  different  tissues  used).  This  will  be  reported  at  the  2003 
International  Ultrasonics  Symposium  in  October. 

Reportable  Outcomes 

Steidl,  C.;  Hut  Yao ;  Phukpattaranont,  P Ebbini,  E.S. ,  "Dual-mode  ultrasound  phased  arrays  for 

noninvasive  surgery:  post-beamforming  image  compounding  algorithms  for  enhanced 


8 


Image-Guided  Surgery  of  Primary  Breast  Cancer  Using  Ultrasound  Phased  Arrays 

Emad  S.  Ebbini.  PI 

visualization  of  therma!  lesions,"  Biomedical  Imaging,  2002.  Proceedings.  2002  IEEE 
International  Symposium  on  ,  7-10  July  2002.  Page(s):  429  -432. 

tesions'"  Proceed ES-'"Nonli^ar  methods  for  visualization  of  HIFU-induced 
luiJ  -  i  a  9s  of  the  2  International  Symposium  on  Therapeutic  Ultrasound  (ISTU2) 

29  July  l  August,  Seattle,  WA,  Editors:  Andrew,  Crum,  and  Vaezy,  2002,  Page(s):  282  -  289. 

dui,Y^ P-l  Ebbinif  ES-'  "Detection  and  mapping  of  thermal  lesions  using 
Volume  2  Ort TiTdonnaJSeDd  Proceedin9s  of  the  2002  IEEE  Ultrasonics  Symposium. 

mFuInL^Kr3"0^'  '#***’.  £s-"N°nlinear  imaging  methods  for  characterization  of 
HIFU-induced  lesions,  Proceedings  of  the  SPIE  vol.  4954,  Thermal  Treatment  of  Tissue:  Energy 
Delivery  and  Assessment  II,  Editor:  T.P.  Ryan,  2003,  Page(s):  183  -  191.  9V 

E  Ebbini'  "Nonlinear  pulse-echo  imaging  methods  for  HIFU-induced  lesion  visualization," 
°f  thS  AC0UStlcal  Society  of  America,  JASA,  vol.  113,  no.  4,  2003,  Page:  2309 

S'  P^a^n^odo  HIFU  arrays,"  3rd  International  Symposium  on  Therapeutic 

Ultrasound  (ISTU3),  Lyon,  France  22  -  25  June,  2003,  manuscript  in  preparation. 

EmadS.  Ebbini  and  Pornchai  Phukpattaranont,  "Ultrasound  imaging  system  and  method  using 
non-linear  post-beamforming  filter,"  Patent  application  to  the  US  Patent  and  Trademark  Office 
(Initial  application  10  May,  2002;  Continuation-in-part  9  May  2003). 

Pornchai  Phukpattaranont  and  Emad  S.  Ebbini,  "Post-beamforming  second-order  Volterra  filter 

imaging",  IEEE  Trans,  on  Ultrasonics,  Ferroelectrics  and  Frequency 


Conclusions 

Our  efforts  during  the  second  year  of  the  grant  were  focused  on  validating  and  characterizing  the 
imaging  capabilities  of  the  64-element  dual-mode  array  prototype.  In  addition,  we  have  generalized 
the  harmonic  processing  of  the  echo  signal  to  quadratic  processing.  This  approach  was  shown  to  be 
superior  in  terms  of  enhancing  the  spatial  and  contrast  resolution  of  lesion  images.  Quadratic  imaging 
was  also  shown  to  produce  results  comparable  to  pulse  inversion  imaging,  an  imaging  method 
designed  to  detect  microbubble  activity  (contrast-agent  imaging).  Compared  to  pulse  inversion 
quadratic  imaging  is  not  sensitive  to  motion  and  has  wider  dynamic  range. 

Our  characterization  of  the  64-element  dual-mode  array  prototype  is  as  follows: 

Therapeutic  performance:  Free  field  focal  intensities  up  to  3500  W/cm2  in  a  region  approximately  4 
cm  m  diameter  around  the  geometric  focus  (with  electronic  scanning). 

Imaging  Performance:  Axial  resolution  of  3  mm,  lateral  resolution  of  1.1  mm.  This  is  opposite  to 
typical  diagnostic  systems,  which  have  better  axial  resolution  (<  1  mm)  than  axial  resolution  (-2  mm). 

ns  is  due  to  the  relatively  narrow  bandwidth  (leading  to  reduced  axial  resolution)  and  the  lar^c 
apeiture  (leading  to  excellent  lateral  resolution)  of  the  dual-mode  array  prototype.  Furthermore,  wc 
nive  show'n  that  the  dynamic  range  of  50  dB  in  an  elliptical  region  around  the  geometric  center  (major 
diameter  6  cm  and  minor  diameter  of  5  cm).  High  contrast  targets  like  large  blood  vessels  are  always 
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reflector SPCCkle  C°ntraSt  targets  smal,er  than  4  mm  are  indistinguishable  from  a  point  specular 

What  does  this  mean?  Our  dual  mode  approach  is  a  viable  approach  for  image  guidance  as  wc  can 
recognize  structures  and  markers  in  the  target  region.  We  have  also  established  unequivocally  that 
lesions  can  be  mapped  accurately  based  on  changes  in  their  echogenecity  during  and  immediately  after 
omiation  This  remains  one  of  the  most  attractive  features  of  our  approach,  especially  when  sin^lc- 
transmit  focusing  is  used.  This  mode  is  unique  because  it  allows  us  to  inspect  the  changes  in  the 
echoes  precisely  at  the  expected  lesion  location  (since  the  imaging  transmit  beam  is  obtained  from  the 
therapeutic  beam  itself).  It  also  allows  for  the  detection  of  any  critical  structures  in  the  path  of  the 
therapeutic  beam.  Therefore,  imaging  with  the  therapeutic  beam  before  lesion  formation  is  probably 
the  most  powerful  tool  for  avoiding  collateral  damage  upon  targeting  specific  tissue  for  damage. 
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DUAL-MODE  ULTRASOUND  PHASED  ARRAYS  FOR  NONINVASIVE  SURGERY* 
POST-BEAMFORMING  IMAGE  COMPOUNDING  ALGORITHMS  FOR  ENHANCED 
VISUALIZATION  OF  THERMAL  LESIONS 
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Absiraci-  A  past-beam fo raiing  nonlinear  compounding  algo¬ 
rithm  for  ultrasonic  imaging  is  presented.  Fundamental  and 
harmonic  image  components  from  beamfromed  radio  frequency 
(RF)  data  are  extracted,  envelope  detected  and  compounded 
using  a  spatial  compounding  functions  (SCFs)  derived  from 
the  transm  it/receive  beamforming  topology  used  in  obtaining 
the  RF  data.  This  is  specially  useful  for  applications  where 
single-transmit  focus  (STF)  imaging  is  used.  In  this  paper,  we 
present  results  from  STF  imaging  experiments  using  a  novel 
dual-mode  phased  array  system  for  image-guided  surgery.  In 
particular,  we  address  the  enhancement  of  the  echogenicity  of 
thermal  lesions  formed  in  ex  vivo  tissue.  It  is  shown  that  new 
nonlinear  image  compounding  algorithm  produces  25  -  30  dB 
enhancement  in  lesion  echogenicity  without  loss  in  spatial  res¬ 
olution.  This  is  to  be  compared  with  a  typical  enhancement 
of  5  dB  achieved  by  standard  echographic  imaging  and  20  dB 
achieved  by  second  harmonic  imaging  alone.  In  addition,  im¬ 
ages  resulting  from  the  new  algorithm  are  virtually  free  of 
beamforming  artifacts  that  can  severely  degrade  the  perfor¬ 
mance  of  2nd  harmonic  imaging. 

i.  INTRODUCTION 

Image  guidance  has  long  been  recognized  as  the  “enabling 
technology”  for  noninvasive  thermal  surgery.  Highly  re¬ 
fined  energy  application  devices  have  been  developed  and 
for  years.  However,  without  reliable  imaging  techniques  for 
visualization  of  thermal  lesions,  noninvasivethennal  surgeiy 
has  failed  to  find  widespread  acceptance  in  the  clinic.  Re¬ 
cently,  image  guidance  methods  based  on  well  established 
imaging  modalities  like  MRJ[3],  CTf5],  or  ultrasound  [4,  8} 
have  been  proposed.  Other  imaging  modalities  are  also  be¬ 
ing  developed  and  may  become  available  in  the  foreseeable 
future[6, 7]. 

An  area  unique  to  ultrasound  that  could  revolutionize 
the  field  of  image-guided  surgery  is  the  development  of  a 
new  generation  of  dual-mode  high-power  phased  array  sys- 

Fundcd  by  Gram  DAMD  17-01-1-330  from  the  US  Army  Medical 
Research  and  Materiel  Command. 


terns  capable  of  both  imaging  and  therapy  [8.  9].  These 
piezocomposite  transducers  can  produce  focal  intensity  lev¬ 
els  needed  for  ablative  and  coagulative  thermal  surgery  with 
high  precision.  Furthermore,  the  operating  bandwidth  of 
such  transducers  allows  for  imaging  the  treatment  region 
with  adequate  image  quality  to  delineate  important  land¬ 
marks  within  and  around  the  target  volume.  With  these  ca¬ 
pabilities,  it  is  possible  to  operate  these  arrays  in  a  “self- 
registration”  mode  whereby  the  imaging  capabilities  of  the 
array  are  utilized  in  characterization  of  the  tissue  response 
precisely  at  the  expected  lesion  location.  This  is  due  to  the 
fact  that  the  beamforming  is  common  to  both  the  imag¬ 
ing  and  therapy  modes.  The  main  challenge  to  this  ap¬ 
proach  is  the  low-contrast  nature  of  speckle-ridden  ultra¬ 
sound  images.  This  paper  addresses  a  new  post  beamform- 
ing  filter  bank  image  reconstruction  algorithm  with  nonlin¬ 
ear  spatiaJy  weighted  compounding  algorithm  to  mitigate 
this  problem.  The  design  of  this  algorithm  is  motivated  by 
the  nonlinear  nature  of  ultrasound  propagation  in  tissue  and 
the  fact  that  thermal  lesions  are  known  to  have  increased 
level  of  harmonic  generation.  Experimental  results  demon¬ 
strate  the  potential  of  the  new  algorithm  in  both  enhancing 
the  lesion  contrast  and  size/shape  analysis. 


2.  IMAGE  FORMATION  MODEL 

Figure  1  summarizes  the  image  acquisition  and  image  for¬ 
mation  model.  A  64-element  array  optimized  for  maximum 
energy  delivery  at  1  MHz  operating  frequency  is  used  for  le¬ 
sion  formation  in  sample  tissue.  Lesions  are  formed  by  fo¬ 
cusing  the  array  at  a  point  within  the  target  and  maintaining 
high-power  output  for  time  intervals  on  the  order  of  seconds 
(1-5  seconds  typical).  The  power  is  interrupted  for  short 
intervals  (milliseconds)  to  acquire  image  data  by  transmit¬ 
ting  short  (jus)  pulses  from  all  64-elements  and  receiving  on 
selected  elements  using  a  matrix  switch.  Once  the  image 
data  set  is  collected,  RF  beamforming  is  performed  to  form 
standard  echographic  images  of  the  target  region. 
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Fig.  I.  Image  Acquisition  and  Image  Reconstruction  Algorithm 


Fig.  2.  Coordinate  system  used  in  the  synthetic  aperture 
imaging  system. 

2.1.  Synthetic  aperture  imaging 

Images  are  obtained  using  a  full  synthetic  aperture  tech¬ 
nique  [10].  The  image  pixel  at  coordinates  (xp,zp)  is  there¬ 
fore  computed  by  (Fig.  2): 

64  64 

/(rP!  *,)  =  ££  At  ■  Bj  ■  SiJ  [(ii,,  +  Rjp)/C] ,  (1) 

i=l  J=1 

where  i  is  the  transmit  element  index,  j  is  the  receive  el¬ 
ement  index.  A*  is  the  transmit  apodization  weight  at  ele¬ 
ment  x,  Bj  is  the  receive  apodization  weight  at  element  j, 
RtP  and  Rjp  are  the  distances  from  the  transmit  and  receive 
elements,  respectively,  from  the  image  pixel,  c  is  the  speed 
of  sound  in  the  medium  being  imaged,  and  st- ,-(t)  is  the 
echo  acquired  when  transmitting  with  element  i  and  receiv¬ 
ing  with  element  j. 

2.2.  Single-TVansmit  Focus  (STF)  Imaging 

The  synthetic  aperture  imaging  algorithm  described  above 
can  be  used  to  produce  the  highest  quality  conventional  im¬ 
ages  that  can  be  expected  from  a  given  array.  However,  due 
to  the  fact  that  transmit  beams  are  synthesized  by  super¬ 
position  of  single-element  transmit  patterns,  the  nonlinear 
interactions  of  the  real-time  transmit  beams  cannot  be  ac¬ 
counted  for  using  this  imaging  algorithm.  Therefore,  we 


have  modified  our  64-channe!  phased  array  driver  to  allow 
for  pulsed  transmission  on  all  64  channels  simultaneously. 
This  allowed  us  to  use  the  full  power  of  the  transmit  beams 
and,  therefore,  observe  their  nonlinear  interactions  with  the 
tissue  media.  The  image  formation  process  due  to  a  STF 
beam  is  a  modified  version  of  Equation  (1)  as  follows: 

64 

I(xp,  zp)  =  ^2  Bj  •  sj  [(Bq  +  Rjj})  /c\ ,  (2) 

where  Sj  (t )  is  the  received  waveform  at  element  j  due  to  the 
transmitted  beam  and  /2<j  is  a  fixed  distance  determined  by 
the  focal  depth  of  the  transmit  beam.  All  other  quantities  in 
Equation  (2)  are  the  same  as  their  counterparts  in  Equation 
( 1 ).  It  should  be  noted  that  beamforming  is  performed  in  the 
RF  domain  using  true  delays,  i.e.,  no  baseband  conversion 
is  done.  This  retains  all  the  frequency  components  in  the 
beamformed  data  for  postbeamfonning  filtering  operations 
described  below. 

23.  Post  Beamforming  Filtering 
Since  beamforming  is  performed  in  the  RF  domain,  all  the 
frequency  components  in  the  received  echo  signals  are  re¬ 
tained.  This  is  important  since  the  echo  signals  can  be  ex¬ 
pected  to  contain  a  mix  of  harmonics  (and  possibly  subhar¬ 
monics)  depending  on  the  nonlinearity  of  the  tissue  medium 
being  imaged.  This  is  especially  true  for  the  STF  imaging 
mode  where  the  transmit  power  is  sufficiently  high  to  pro¬ 
duce  harmonic  components  in  nonlinear  tissue  media.  Post- 
beamforming  filtering  can  be  used  to  isolate  specific  har¬ 
monic  components  and/or  enhance  axial  resolution  if  the 
SNR  of  the  system  is  sufficiently  high.  Algorithms  for  post¬ 
beamfonning  fiherbank  image  reconstruction  for  pulse-echo 
ultrasound  can  be  found  in  [II].  For  the  purposes  of  this 
paper,  the  filter  bank  is  formed  of  bandpass  filters  with  cen¬ 
ter  frequencies  at  a  set  of  preselected  harmonics  (or  subhar¬ 
monics). 

2.4.  Nonlinear  Compounding 

Ultrasound  propagation  in  tissue  media  is  inherently  nonlin¬ 
ear.  RF  echo  signals  obtained  using  modem  scanners  with 
wideband  transducers  carry  harmonic  components  that  are 
generated  by  the  medium  nonlinearity,  i.e.,  they  are  not  pan 
of  the  originally  transmitted  imaging  pulses.  Since  the  tis¬ 
sue  nonlinearity  parameter  exhibits  higher  contrast  between 
vanous  tissue  compared  to  tissue  reflectivity,  imaging  this 
parameter  can  lead  to  higher  contrast  images  of  soft  tissue. 
For  a  variety  of  physical  reasons,  it  is  well  known  that  ther¬ 
mal  lesions  generate  higher  harmonics  at  levels  higher  than 
normal  tissues.  By  careful  design  of  the  transmit/receive 
beamforming  and  the  post  beamforming  reconstructions  fil¬ 
ters,  high  contrast  images  of  lesions  can  be  formed  from 
ultrasound  images  formed  at  the  fundamental  and  some  of 
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its  higher  harmonics.  In  general,  an  optima]  compound¬ 
ing  rule  must  be  derived  based  on  statistical  model  of  the 
imaged  region  and  the  expected  lesion  size/location.  Fur¬ 
thermore,  the  quality  of  the  beamforming  at  different  har¬ 
monics  or  scales  will  be  spatially  heterogeneous.  For  exam¬ 
ple,  since  the  array  sampling  function  (in  wavelengths)  is 
coarser  at  the  higher  harmonics,  images  formed  at  these  fre¬ 
quency  will  have  high  contrast  along  the  axis  of  the  trans¬ 
mit  beam,  but  the  contrast  deteriorates  quickly  away  from 
the  mam  axis.  Furthermore,  since  the  harmonics  are  orders 
of  magnitudes  smaller  than  the  fundamental,  compounding 
is  performed  with  log-compressed  images  after  appropriate 
scaling.  Therefore,  with  reference  to  Figure  1,  the  com¬ 
pounding  is  performed  in  two  stages: 

1.  The  outputs  of  the  filterbank  are  envelope  detected, 
scaled,  and  log  compressed.  The  scaling  of  the  indi¬ 
vidual  harmonics  can  be  derived  from,  for  example, 
the  specificity  or  contrast-to-tissueratio  (CTR): 

CTfl=10,OE'"(iraf)  °> 

where  ||yc,  ((2  and  Hyr;  H2  are  the  I2  norms  of  the  tth 
harmonic  components  from  the  contrast  and  normal 
tissue  regions,  respectively.  These  regions  arc  easily 
identified  under  various  imaging  conditions.  For  in¬ 
stance,  for  the  application  described  in  this  paper,  the 
contrast  region  is  the  expected  location  of  the  ther¬ 
mal  lesion  (often  visible  on  the  standard  echographic 
image). 

2.  A  spatially-weighted  sum  of  the  harmonic  compo¬ 
nents  is  performed.  The  weighting  function  is  derived 
from  the  spatial  contrast  function  (SCF)  of  the  imag¬ 
ing  array  at  the  different  harmonics. 

Assuming  that  the  envelope -detected,  normalized  and  log- 
compressed  images  obtained  from  the  filterbank  are  given 
by  A  ?  ^2:  •  •  * » //v.  then  the  compounded  image  is  given  by: 

AT 

I(x>y)  =  (4) 

issZ 

where  is  a  weighting  function  reflecting  the  relative  en¬ 
ergy  at  the  harmonic  with  respect  to  the  fundamental  at  the 
(expected)  lesion  location  and  5i(x,y)  is  the  SCF  associ¬ 
ated  with  the  ith  harmonic.  This  is  can  be  obtained  based 
on  simulations  of  pulsed  or  CW  calculations  of  the  follow¬ 
ing  ratio: 


(x  y)  =  Lik  Br'  ^  yoi^)dx 


where  the  subscripts  sl  and  \fi  refer  to  sidelobe  and  main- 
lobe,  respectively,  and  £r.  (ar?  y;  s')  is  the  ith  harmonic  (dy¬ 
namic)  receive  focus  at  (x,  y)  and  Bt(x,  y0;  x')  is  the  trans¬ 
mit  focus  pattern  (at  the  fundamental).  Note  that  the  SCF 


is  a  measure  of  the  quality  of  the  transmit/receive  partem  at 
the  ith  harmonic.  The  higher  the  SCF  the  more  confidence 
we  have  in  the  received  data  at  the  ith  harmonic. 

3.  RESULTS 

Figure  3  gives  a  typical  result  of  a  visualization  experiment 
using  the  dual-mode  array  and  the  image  compounding  al¬ 
gorithm  described  above.  The  figure  shows  images  obtained 
before  and  after  lesion  formation  in  an  ex- vivo  liver  tissue 
samples.  The  expected  lesions  are  cigar-shaped  with  length 
of  10  mm  and  width  of  2  -  3  mm.  All  images  are  shown 
at  40  dB  dynamic  range  and  have  the  same  axial  and  lateral 
extent  The  lesion  was  formed  at  the  geometric  center  of  the 
array  and  was  expected  between  90  and  100  mm  in  the  axial 
direction.  All  images  show  the  front  surface  of  the  sample  at 
80  mm  and  the  sponge  backing  phantom  at  120  mm.  Three 
pairs  of  images  are  shown,  one  before  (left)  and  one  after 
the  lesion  formation.  The  images  on  top  of  each  figure  are 
formed  at  the  fundamental  ( 1  MHz)  while  those  at  the  center 
of  each  figure  are  formed  at  the  2nd  harmonic  (2  MHz).  The 
pair  at  the  bottom  shows  tbe  compound  image  using  Equa¬ 
tion  4  and  the  SCFs  shown  in  Figure  4.  Each  pair  of  images 
is  displayed  on  a  log  scale  and  normalized  such  that  the  0 
dB  corresponds  to  the  maximum  of  the  image  on  the  RHS, 
i.e.,  the  image  of  the  target  after  lesion  formation.  While 
both  the  fundamental  and  harmonic  images  show  enhanced 
lesion  contrast,  the  harmonic  images  show  a  net  increase  in 
contrast  by  22  dB  along  the  direction  of  the  transmit  beam. 
On  the  other  hand,  the  net  increase  in  contrast  at  the  funda¬ 
mental  is  about  7  dB.  In  addition,  the  spatial  definition  of 
the  lesion  in  the  harmonic  images  is  superior  to  the  funda¬ 
mental.  However,  while  the  harmonic  images  are  superior 
both  in  terms  of  definition  and  contrast  near  the  main  axis  of 
the  transmit  beam,  the  fundamental  images  better  represent 
the  tissue  state  off-axis.  The  compound  images  show  a  net 
increase  in  the  lesion  contrast  on  the  order  of  35  dB  both 
axially  and  laterally.  This  result  can  be  understood  in  light 
of  a  close  examination  of  the  SCR s  shown  in  Figure  4.  It  is 
also  interesting  to  note,  however,  that  the  compound  image 
still  shows  major  off-axis  objects,  e.g.,  blood  vessel  at  axial 
distance  of  1 10  mm  and  lateral  distance  of  5  rum.  That  is, 
tbe  compounding  algorithm  removes  the  beamforming  side- 
lobe  artifacts,  but  not  major  scartcrers  in  the  vicinity  of  the 
transmit  beam.  This  is  one  of  the  main  objectives  of  STF 
imaging  described  in  this  paper. 

4.  CONCLUSIONS 

We  have  shown  that  post  beamforming  filterbank  recon¬ 
struction  of  ultrasound  images  ai  selected  frequencies  sen¬ 
sitive  to  harmonic  generation  in  nonlinear  media  produces 
images  appropriate  for  compounding.  The  results  of  apply¬ 
ing  the  nonlinear  compounding  algorithm  also  demonstrate 
a  significant  increase  in  lesion  contrast  over  enhancement 
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Fig.  3.  Images  before  (left)  and  after  lesion  formation.  AU  images 
are  shown  at  40  dB  dynamic  range. 


achieved  at  either  the  fundamental  or  second  harmonic.  It  is 
important  to  state  that  this  contrast  enhancement  is  achieved 
without  loss  in  spatial  resolution.  This  due  to  the  fact  that 
the  bandpass  filters  used  in  separating  the  fundamental  and 
harmonic  components  have  wider  bandwidth  than  the  data 
bandwidth  at  these  harmonics.  The  results  shown  in  this  pa¬ 
per  suggest  that  dual-mode  arrays  with  high  contrast  imag¬ 
ing  capability  suitable  for  real-time  thermal  lesion  visual¬ 
ization  are  feasible.  This  leads  to  a  most  powerful  paradigm 
for  image  guided  surgery  where  the  therapeutic  and  image- 
guidance  coordinate  systems  are  inherently  registered. 
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jU^r^Uj^uta^  imaging  methods  (like  pulse  inversion  [1]  and  quadratic  imagine 

titsS  R°dhfr  Volle[™  fi|ter  f-J)  are  used  in  visualization  of  lesion  fonnation  in  frcshlv  exetscd 
tissue.  Both  of  these  methods  arc  more  sensitive  to  nonlinear  echoes  (e.g..  due  to  micro-bubbles) 
han  standard  B-mode  imaging.  While  all  three  methods  typically  show  increased  echogenicity  at 
the  lesion  location,  the  nonlinear  methods  exhibit  more  localized  echo  enhancement  than  B-mode 
imaging.  Therefore,  nonlinear  methods  are  better  suited  to  lesion  mapping  for  purposes  of  image 
guidance.  Quadratic  images  have  the  added  advantage  of  a  significant  increase  in  image  dynamic- 
range  and  noise  reduction  (a  major  limitation  of  pulse  inversion  imaging).  The  results  shown  in 
this  report  continue  to  support  the  hypothesis  that  micro-bubbles  play  an  important  role  of  lesion 
formation  Furthermore,  the  presence  of  microbubbles  provides  significant  opportunity  for  mapping 
the  treated  tissue  and  potentially  characterizing  the  nature  of  damage. 


INTRODUCTION 

We  have  used  a  dual-mode  array  described  in  |3]  to  form  HIFU-induced  thermal  lesions 
in  freshly  excised  degassed  tissue  under  a  variety  of  normal  exposure  and  over  exposure 
conditions.  Single-transmit  focus  images  were  collected  for  over  100  lesions  before  and 
alter  lesion  formation.  These  images  have  consistently  shown  5  -  7  dB  enhancement  in 
the  echogenicity  from  the  lesion  location  in  the  standard  echographic  images.  These  re¬ 
sults  were  much  more  consistent  than  the  reported  “flashes"  on  the  B-scarfimages  when 
diagnostic  ultrasound  systems  are  used  to  monitor  HIFU  lesion  formation.  Motivated  by 
the  excellent  investigation  by  P.  P.  Lele  reported  in  [4J,  we  hypothesized  that  this  change 
in  echogenicity  is  due  to  stable  microbubbles  that  can  occur  even  at  low  insonation  lev¬ 
els.  Lele  found  that  subharmonic  emission  due  to  microbubbles  showed  a  monotonic 
increase  with  intensity  from  150  mW/cm2  to  1500  W/cm2  without  a  distinct  threshold 
for  emission  (measurements  done  in  vitro  and  in  vivo  at  2.7  and  1.8  MHz).  The  consis¬ 
tency  of  the  increase  in  echogenicity  at  the  lesion  may  be  explained  by  the  fact  that  the 
microbubbles  may  already  be  resonant  at  the  imaging  frequency  (same  as  the  therapeutic 
HIFU  beam  when  the  dual-mode  array  is  is  used),  perhaps  a  result  of  rectified  diffusion. 

The  standard  echographic  images  at  the  fundamental,  however,  offer  limited  contrast 
enhancement  due  to  the  speckle  phenomenon.  Therefore,  they  could  not  provide  a 
reliable  method  for  mapping  the  boundaries  of  HIFU-induced  lesions.  This  lead  us  to 
try  to  exploit  the  nonlinear  nature  of  the  microbubbles  to  enhance  the  visualization  and 
mapping  of  thermal  lesions.  The  idea  was  that,  if  microbubbles  are  indeed  present  at 
the  lesion  location,  they  will  generate  nonlinear  echoes  that  may  be  better  suited  for 


mapping.  We  have  initially  investigated  second  harmonic  (SH)  imaging  as  a  means  of 
enhancing  the  lesion  contrast  for  improving  the  visualization  of  these  images.  SH  imaizes 
or  thermal  lesions  have  shown  increase  in  the  contrast  on  the  order  of  22  -  25  dB.  but  with 
decreased  dynamic  range  of  the  resulting  images  [3J.  A  post-beamforming  nonlinear 
compounding  algorithm  was  shown  to  improve  the  contrast  in  lesion  echogenicity  to 
~r  '  35.dB  WI*out  loss  dynamic  range  [5J.  This  was  achieved  by  compoundin'’ 
the  fundamental  and  the  SH  images  using  spatial  compounding  functions  based  on 
the  receive  beamforming  characteristics  of  the  dual-mode  array  at  the  fundamental 
311  SH  frecluencies-  In  this  paper,  we  use  a  commercial  imaging  scanner  with 
modifications  to  allow  pulse  inversion  imaging  in  .addition  to  standard  B-mode  imagine 

[°lrrihe]  visualizatl0n  of  freshly  excised  tissue  before,  during,  and  after  the  formation  of 
HIFU  lesions. 


NONLINEAR  IMAGING  METHODS 
Pulse  Inversion  Imaging 

This  method  was  recently  introduced  by  Bums  and  coworkers  [6J  for  enhancing 
contrast  echoes  in  contrast-assisted  pulse-echo  imaging.  The  basic  implementation  of 
this  method  entails  the  use  of  two  pulses  per  image  line.  These  pulses  are  carefully 
designed  such  that,  p2(t)  =  -p,(r  -  7J,  where  TL  is  some  appropriate  delay  (on  the 
order  of  the  pulse-echo  time  from  the  maximum  depth  of  interest).  Summing  the  echoes 
resulting  from  the  two  pulses  eliminates  the  odd-harmonic  components  from  the  echo 
signal  (including  the  fundamental)  while  doubling  the  even-harmonic  (mostly  second) 
components.  This  method  currently  represents  the  leading  approach  for  contrast-agent 
imaging,  especially  with  low  concentration  and/or  very  low  transmit  signals  to  minimize 
the  generation  of  tissue  nonlinearity. 


Quadratic  Imaging  Based  On  SVF 

We  have  recently  developed  a  new  nonlinear  imaging  system  based  on  the  SVF  [2J. 
This  has  a  number  of  advantages  when  compared  to  PI  imaging  (e.g.  requires  only  a 
single  pulse  per  line  and  increased  dynamic  range). 


Second-Order  Volterra  Model 

Results  from  [2]  have  shown  the  validity  of  a  second-order  Volterra  filer  as  a  model 
for  pulse-echo  ultrasound  imaging  data  from  tissue  mimicking  media.  In  this  section,  the 
decomposition  of  received  echo,  i.e.,  output  sequences  only,  into  linear  and  quadratic 
components  by  using  least-squares  approach  of  second-order  Volterra  model  will  be 
considered  and  the  detail  of  algorithm  implementation  to  pulse-echo  ultrasound  imaaimt 
will  be  stated. 


Signal  Separation  Model 

The  algorithm  described  in  this  section  is  adapted  from  [7J.  The  response  of  a  quadrat- 
ically  nonlinear  system,  y(n  +  1).  can  be  predicted  by  a  second-order  Volterra  model  of 
past  m  values  as  follows: 

y(n  +  1 )  =  yL{n+  l)+ye(n+  1) 

=  v(«  - +  £  ^y(n-y')y(«-A-)Ae(y./:)  +  e(n).  (1) 

*=o  y=0  k= j 

where  hL{i)  is  linear  filter  coefficients.  hQ{j,k )  represents  quadratic  filter  coefficients 
and  e(n)  is  a  modelling  error  and/or  a  measurement  noise  which  is  assumed  to  be  an 
independent,  identically  distributed(i.i.d)  random  variable  with  zero  mean.  That  is,  if 
the  model  coefficients  are  known,  the  echo  signal  can  be  decomposed  into  linear  and 
quadratic  components.  The  latter  can  be  expected  to  better  represent  the  quadratic  re¬ 
sponse  ol  the  system  than,  say,  the  second  harmonic  component.  The  model  coefficients 
can  be  obtained  by  setting  the  linear  and  quadratic  prediction  problem  in  Equation  1 
to  form  a  set  of  linear  equations.  Recognizing  that  the  output  is  linear  in  terms  of  the 
(unknown)  model  coefficients,  one  obtains  a  matrix  equation  of  the  form: 

f=Gh  +  e.  (2) 

where  the  vector  f,  the  matrix  G  and  the  error  vector  e  are 

f  =  [y(«+l),.v(n  +  2) . y(«  +  Z.)]r 

G=  [y(«)-y(«+ 1) — y(rt+T- i)]r 

e  =  [£(«),£(«+  1), ...,e{n  +  L-  l)]r. 
where  the  data  vector,  y,  is  given  by: 


y(«)  =  [y(«).y(n-  l),y(/i-2),...,y(/i-m-l-  1). 
y2(n),y(n)y(n-  1 ),...,  y2(n-m  +  l)]r 

and  the  filter  coefficient  vector,  h,  is  given  by: 

h  =  [hL{0).hL(\).hL(2) . hL(m  -  1), 

/Jq(0,0),/z^(0,  1 hg^m  -  l.m  -  1)]T. 

The  details  of  the  solution  for  the  coefficients  of  the  SVF  model  can  be  found  in 
[8].  Briefly,  a  minimum-norm  least-squares  solution  of  (2)  is  obtained  using  truncated 
singular  value  decomposition,  TSVD.  To  assess  the  performance  of  the  signal  separation 
model  in  enhancing  the  lesion  visualization,  we  compute  the  contrast-to-tissue  ratio: 

Hygcj\ 
llygrlli  / 


CTR  =  1 0  log,0 


(3) 


where  ||y<v7||2  and  ||yCT||2  are  the  l2  norms  of  the  quadratic  components  from  the  lesion 
and  normal  tissue  regions,  respectively.  These  regions  are  easily  identified  under  various 
imaging  condit.ons.  For  instance,  for  the  application  described  in  this  paper,  the  contrast 

region  is  the  expected  location  of  the  thermal  lesion  (often  visible  on  the  standard 
echographic  image). 


Implementation 


The  coefficients  of  the  linear  and  quadratic  components  of  the  SVF  model 
from  the  beamformed  RF  data.  The  implementation  steps  are  as  follows: 


are  obtained 


1 .  Using  the  standard  echographic  image,  select  a  beamformed  RF  data  segment  from 
the  expected  lesion  location. 

2.  Form  the  linear  systems  of  equations  according  to  (2). 

3.  Define  contrast  region  (within  the  lesion)  and  normal  tissue  region  for  the  compu¬ 
tation  of  the  mean-square  error  MSE  and  CTR.  r 

4.  Solve  systems  of  linear  equations  by  using  TSVD  regularization  method. 

5.  Apply  second-order  Volterra  filter  to  the  beamformed  RF  data  throughout  the  pulse- 
echo  ultrasound  image. 

6.  The  quadratic  component  from  the  SVF  can  be  displayed  as  a  separate  image  or 
appropriately  compounded  with  the  linear  component. 


EXPERIMENTAL  SETUP 

Figure  1  shows  a  simple  arrangement  for  the  formation  of  HIFU  lesions  in  freshly 
excised  and  degassed  porcine  livers  samples.  The  therapy  transducer  is  a  1.5  MHz 
single-element  spherical-shell  transducer  with  a  radius  of  curvature  equal  to  its  diameter 
and  equal  to  63.5  mm  (Etalon.  Lizton,  Indiana).  The  transducer  is  fixed  to  the  back  of 
a  small  tank  as  shown  and  driven  by  a  power  amplifier  (ENI,  Rochester,  NY)  and  a 
programmable  function  generator.  This  assembly  can  be  used  in  generating  a  variety  of 
amplitude-modulated  HIFU  bursts  from  tens  of  milliseconds  to  several  seconds  lon«  and 
intensities  up  to  3000  W/cm2  (conservative  estimate). 

Real-time  imaging  is  performed  using  a  modified  Technos  MP  system  from  ESAOTE. 
enoa,  Italy.  The  system  is  modified  to  allow  imaging  in  pulse  inversion  mode  in 
addition  to  normal  B-mode  imaging.  In  addition,  a  hardware  module  for  capturing  high- 
quality  beamformed  RF  data  allows  us  to  capture  and  upload  up  to  60  seconds  °of  full 
rame  data  with  a  specified  frame  rate.  A  CA  421  convex  probe  was  used  in  acquiring 
image  data  for  this  paper.  Image  data  was  acquired  in  pulse  inversion  (PI)  mode  with 
a  --cycle  transmit  pulse  centered  at  1.57  MHz.  The  imaging  transducer  was  aligned  so 

that  the  image  plane  is  x  -  z  plane  (to  allow  imaging  the  lesion  along  the  axis  of  the 
therapy  transducer. 


Therapy  Transducer  (Thx  Tx) 


Imaging  Transducer  (lx  Tx) 


Acoustic  Window 


FIGURE  I  Experiment  setup  used  lor  formation  of  HIFU  lesions  using  a  single-clement  transducer 

attached  to  the  of  the  tank  with  an  imaging  probe  monitoring  a  cross  section  of  the  lesion  through  an 
acoustic  window.  c 


RESULTS  AND  DISCUSSION 

The  experimental  setup  shown  in  Figure  l  was  used  in  obtaining  images  of  ex  vivo  tissue 
samples  before,  during  and  after  thermal  lesion  formation  with  the  spherical  therapy 
transducer.  This  setup  was  designed  to  allow  comparisons  with  our  a  dual-mode  64- 
element  ultrasound  phased  array  operating  at  I  MHz  described  in  [3J.  Results  shown  in 
[3J  confirmed  that  echoes  from  the  lesion  location  exhibited  increased  levels  of  second 
harmonic  generation  even  at  normal  exposure  conditions.  In  this  paper,  we  show  imaging 
results^  from  a  single-shot  slightly  overexposure  condition  of  a  5  second  pulse  at  1800 
W/cm  .  This  exposure  typically  produces  a  cigar-shaped  lesion  with  average  length  of 
10  mm  and  slight  widening  of  the  lateral  extent  of  the  lesion  at  the  base  (nearer  to  the 
therapeutic  transducer).  This  is  part  of  a  study  for  characterization  of  different  imaging 
modes  at  a  full  range  of  exposure  conditions. 

Figure  2  shows  the  RF  data  along  a  line  through  the  lesion  before  lesion  formation 
(upper  left).  The  figure  also  shows  spectrograms  of  the  RF  data  (showing  the  frequency 
content  of  the  RF  echoes,  lower  left).  Echoes  from  the  tissue  sample  begin  at  depth  135 
mm.  The  results  show  that  the  echoes  before  lesion  formation  cover  the  bandwidth  of  the 
CA4121  probe  (fundamental  at  1.57  MHz  and  some  2nd  harmonic).  This  is  confirmed 
by  the  corresponding  PI  data  shown  on  the  right  hand  side,  which  shows  the  second 
harmonic  data  just  above  the  noise  level. 

On  the  other  hand,  the  spectrogram  ol  the  RF  data  shown  in  Figure  3  from  the  same 
direction  after  lesion  formation  show  strong  fundamental  and  2nd  harmonic  components 
at  the  lesion  location  (145  mm).  This  is  also  very  clear  in  the  PI  data  which  shows  a 
strong  2nd  harmonic  at  the  same  location.  In  addition,  the  PI  data  consistently  showed 
a  significant  component  at  2.2  MHz  (believed  to  be  an  ultra  harmonic  component). 
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ABSTRACT 

It  has  long  been  recognized  that  thermal  lesions  formed  us¬ 
ing  high-intensity  focused  ultrasound  (HIFU)  exhibit  non¬ 
linear  behavior  that  can  be  detected  in  pulse-echo  ultrasound. 
Second  harmonic  imaging  of  freshly  formed  thermal  lesions 
have  consistently  shown  significant  enhancement  in  their  vi¬ 
sualization  confirming  this  nonlinear  behavior.  In  this  pa¬ 
per,  we  describe  a  post-beamforming  nonlinear  filtering  al¬ 
gorithm  based  on  second-order  Volterra  filter  (SVF)  model 
that  separates  the  linear  and  quadratic  components  of  the 
echo  signal  leading  to  significant  enhancement  of  lesion  vi¬ 
sualization.  Images  from  ex  vivo  tissue  samples  are  shown 
to  demonstrate  the  level  of  contrast  enhancement  achieved 
with  the  SVF-based  quadratic  filter  compared  with  standard 
echo  and  2nd  harmonic  imaging  results. 

I.  INTRODUCTION 

Since  early  2001 ,  we  have  used  a  dual-mode  array  described 
m  [1]  to  form  HIFU-induced  thermal  lesions  in  freshly  ex¬ 
cised  degassed  tissue  under  a  variety  of  normal  exposure 
and  over  exposure  conditions.  Singie-transmit  focus  im¬ 
ages  were  collected  for  over  100  lesions  before  and  after 
lesion  formation.  These  images  have  consistently  shown  5 
-  7  dB  enhancement  in  the  echogenicity  from  the  lesion  lo¬ 
cation  in  the  standard  echographic  images.  These  results 
were  much  more  consistent  than  the  reported  “flashes”  on 
the  B-scan  images  when  diagnostic  ultrasound  systems  are 
used  to  monitor  HIFU  lesion  formation.  Motivated  by  the 
excellent  investigation  by  P.  P.  Lelc  reported  in  (2J,  we  hy¬ 
pothesized  that  this  change  in  echogenicity  is  due  to  sta¬ 
ble  microbubbles  that  can  occur  even  at  low  insonation  lev¬ 
els.  Lele  found  that  subharmonic  emission  due  to  microbub- 
bles  showed  a  monotonic  increase  with  intensity  from  150 
mW/cm2  to  1500  W/cm2  without  a  distinct  threshold  for 
emission  (measurements  done  in  vitro  and  in  vivo  at  2.7  and 
1 .8  MHz).  The  consistency  of  the  increase  in  echogenicity 
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at  the  lesion  may  be  explained  by  the  fact  that  the  microbub- 
bles  may  already  be  resonant  at  the  imaging  frequency  (same 
as  the  therapeutic  HIFU  beam  when  the  dual-mode  array  is 
is  used),  perhaps  a  result  of  rectified  diffusion. 

The  standard  echographic  images  at  the  fundamental , 
however,  offer  limited  contrast  enhancement  due  to  the  speckle 
phenomenon.  Therefore,  they  could  not  provide  a  reliable 
method  for  mapping  the  boundaries  of  HIFU-induccd  le¬ 
sions.  This  lead  us  to  try  to  exploit  the  nonlinear  nature 
of  the  microbubbles  to  enhance  the  visualization  and  map¬ 
ping  of  thermal  lesions.  The  idea  was  that,  if  microbuhblcs 
arc  indeed  present  at  the  lesion  location,  they  will  generate 
nonlinear  echoes  that  may  be  better  suited  for  mapping.  We 
have  initially  investigated  second  harmonic  (SH)  imaging 
as  a  means  of  enhancing  the  lesion  contrast  for  improving 
the  visualization  of  these  images.  SH  images  of  thermal  le¬ 
sions  have  shown  increase  in  the  contrast  on  the  order  of 
22  -  25  dB,  but  with  decreased  dynamic  range  of  the  re¬ 
sulting  images  [1],  A  post-beamforming  nonlinear  com¬ 
pounding  algorithm  was  shown  to  improve  the  contrast  in 
lesion  echogenicity  to  30  -  35  dB  without  loss  in  dynamic 
range  [3].  This  was  achieved  by  compounding  the  funda¬ 
mental  and  the  SH  images  using  spatial  compounding  func¬ 
tions  based  on  the  receive  beamforming  characteristics  of 
the  dual-mode  array  at  the  fundamental  and  the  SH  frequen¬ 
cies.  In  this  paper,  we  describe  a  post-beamforming  filter¬ 
ing  algorithm  for  separating  the  linear  and  quadratic  com¬ 
ponents  of  the  beamformed  data  based  on  the  SVF  This 
approach  is  superior  to  all  of  the  algorithms  based  on  SH 
imaging.  It  greatly  enhances  lesion  to  tissue  contrast  with¬ 
out  any  loss  in  spatial  resolution.  In  addition,  it  enhances 
the  dynamic  range  of  the  image  thus  greatly  improving  both 
detecting  and  mapping  of  HIFU-induced  lesions,  even  for 
volumetric  lesions. 


2.  IMAGE  FORMATION 

Figure  1  summarizes  the  image  acquisition  and  image  for¬ 
mation  model.  A  64-elemcnt  array  optimized  for  maximum 
energy  delivery  at  1  MHz  (Imasonic,  Besan<;on,  France) 
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m  values  as  follows: 


Fig.  1.  Imaging  Setup 


is  used  for  lesion  formation  in  sample  tissue.  Lesions  are 
formed  by  focusing  the  array  at  a  point  within  the  target 
and  maintaining  high-power  output  for  time  intervals  on  the 
order  of  seconds  (1-5  seconds  typical).  The  power  is  in¬ 
terrupted  for  short  intervals  (milliseconds)  to  acquire  image 
data  by  transmitting  short  (ps)  pulse  from  all  64-elements 
and  receiving  on  selected  dements  using  a  matrix  switch 
Once  the  image  data  set  is  collected,  RF  beamforming  is 
performed  to  form  standard  cchographic  images  of  the  tar¬ 
get  region.  Figure  1  also  shows  a  general  post-beamforming 
filter  bank  for  analysis  of  beamformed  data.  For  exam¬ 
ple,  the  filter  bank  could  be  designed  to  perform  harmonic 
decomposition  of  beamfromed  data  (e.g.,  fundamental  and 
2nd  harmonic  imaging).  Alternatively,  the  filter  bank  can 
be  designed  to  separate  the  linear  and  quadratic  (nonlin¬ 
ear)  components  based  on  the  2nd  order  Volten-a  model  de¬ 
scribed  below. 


3.  SECOND-ORDER  VOLTERRA  MODEL 

Results  from  [4]  have  shown  the  validity  of  a  second-order 
Volterra  filer  as  a  model  for  pulse-echo  ultrasound  imaging 
data  from  tissue  mimicking  media.  In  this  section,  the  de¬ 
composition  of  received  echo,  i.e.,  output  sequences  only, 
into  linear  and  quadratic  components  by  using  least-squares 
approach  of  second-order  Vblterra  model  will  be  considered 
and  the  detail  of  algorithm  implementation  to  pulse-echo  ul¬ 
trasound  imaging  will  be  stated. 


3.1.  Signal  Separation  Model 

The  algorithm  described  in  this  section  is  adapted  from  f5J. 
The  response  of  a  quadraiically  nonlinear  system,  y(n  -f  1), 
can  be  predicted  by  a  second-order  Volterra  mode!  of  past 


y(n  4-  1}  —  yi{n  +  1)  +  yQ(n  +  1) 

m—  1 

]P  y(n  -  i)hL(i) 

t=0 

m—  1 m  —  1 

+  H  V(n  -  •?>("  -  k)hQ{j,  k)  +  Sin).  1 1 ) 

i=o  k=, 

where  hL(i)  is  linear  filter  coefficients.  hQ{j,  k)  represents 
quadratic  filter  coefficients  and  £{tl )  is  a  modelling  error 
and/ora  measurement  noise  which  is  assumed  to  be  an  inde¬ 
pendent,  identically  distributed(i.i.d)  random  variable  with 
zero  mean.  That  is,  if  the  model  coefficients  are  known,  the 
echo  signal  can  be  decomposed  into  linear  and  quadratic 
components.  The  latter  can  be  expected  to  better  represent 
the  quadratic  response  of  the  system  than,  say,  the  second 
harmonic  component.  The  model  coefficients  can  be  oh- 
tained  by  setting  the  linear  and  quadratic  prediction  problem 
in  Equation  1  to  form  a  set  of  linear  equations.  Recognizing 
that  the  output  is  linear  in  terms  of  the  (unknown)  model 
coefficients,  one  obtains  a  matrix  equation  of  the  form: 

f  =  Gh  +  e,  (2) 

where  the  vector  f,  the  matrix  G  and  the  error  vector  «  are 

f  =  [y(«  +  l),y(n  +  2> . y(n  +  L)  jT 

G=  [y(n),y(n  +  I),...,y(n  +  L-l)lT 

e  =  [e(n),e(n  +  l),...,e(n  +  Z.-  l)]T. 

where  the  data  vector,  y,  is  given  by: 

y(n)  =  [y  W ,  y(n  -  1), y(n  -  2), . . . ,  y (n  ~  m  +  1 ), 

y2M,y(n)y(n  -  1) . y2(n  -  m  4-  1)]T 

and  the  filter  coefficient  vector,  h,  is  given  by: 

b  =  M),  Ml),  M2), ,  hL{m  -  1), 

Aq(0,0),  A<5(0, 1),. . . ,  h.Q(m  -  1,  m  -  1)]T. 

The  details  of  the  solution  for  the  coefficients  of  the  SVF 
model  can  be  found  in  [6].  Briefly,  a  minimum-norm  least- 
squares  solution  of  (2)  is  obtained  using  truncated,  singular 
value  decomposition ,  TSVD.  To  assess  the  performance  of 
the  signal  separation  model  in  enhancing  the  lesion  visual¬ 
ization,  we  compute  the  contrast-to-tissue  ratio: 


CTR  =  I0JogI0 


lly  qc\\1\ 
IlyQT'IIl  / 


(3) 


where  HyQclU  and  IlyQTlh  are  the  I2  norms  of  the  quadratic 
components  from  the  lesion  and  normal  tissue  regions,  re¬ 
spectively.  These  regions  are  easily  identified  under  vari¬ 
ous  imaging  conditions.  For  instance,  for  the  application 
described  in  this  paper,  the  contrast  region  is  the  expected 
location  of  the  thermal  lesion  (often  visible  on  the  standard 
echographic  image). 
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3.2.  Implementation 

The  cocfficiems  of  the  linear  and  quadratic  components  of 
the  SVF  model  are  obtained  from  the  beamformed  RF  data 
t  he  implementation  steps  are  as  follows: 


1  Using  the  standard  echographic  image,  select  a  beam- 

formed  RF  data  segment  from  the  expected  lesion  lo- 
cation. 

2.  Form  the  linear  systems  of  equations  according  to  (2). 

3.  Define  contrast  region  (within  the  lesion)  and  normal 
tissue  region  for  the  computation  of  the  mean-square 
error  MSE  and  CTR, 

4.  Solve  systems  of  linear  equations  by  using  TS  VD  reg¬ 
ularization  method.  B 

5.  Apply  second-order  Volterra  filter  to  the  beamformed 
RF  data  throughout  the  pulse-echo  ultrasound  image. 

6.  The  quadratic  component  from  the  SVF  can  be  dis¬ 
played  as  a  separate  image  or  appropriately  compounded 
with  the  linear  component. 


4.  RESULTS  AND  DISCUSSION 

The  experimental  setup  shown  in  Figure  1  was  used  in  ob¬ 
taining  images  of  ex  vivo  tissue  samples  before  and  after 
thermal  lesion  formation  with  a  dual-mode  64-elemem  ul¬ 
trasound  phased  array  operating  at  1  MHz.  The  standard 
echographic  images  were  formed  by  single  transmit  beams 
[  I  a  ong  the  main  axis  of  the  array  and  dynamic  receive 
focusing  at  each  pixel  in  the  image.  Results  shown  in  [1] 
confirmed  that  echoes  from  the  lesion  location  exhibited  in¬ 
creased  levels  of  second  harmonic  generation  even  at  nor¬ 
mal  exposure  conditions.  Figure  2  shows  the  RF  data  along 
central  axis  of  the  dual-mode  array  before  and  after  lesion 
formation  at  normal  exposure  level  (850  W/cm2  for  4  sec¬ 
onds).  The  figure  also  shows  spectrograms  of  the  RF  data 
(showing  the  frequency  content  of  the  RF  echoes  in  the  ax¬ 
ial  direction).  The  results  show  that  the  echoes  before  le¬ 
sion  formation  are  centered  at  I  MHz  with  no  evidence  of 
2nd  harmonic  component  in  the  tissue  region  (from  80  to 
120  mm).  On  the  other  hand,  the  spectrogram  of  the  RF 
data  from  the  same  direction  after  lesion  formation  shows  a 
strong  2nd  harmonic  component  at  the  lesion  location  (90 
-  100  mm).  Figure  3  shows  the  RF  data  along  central  axis 
of  the  dual-mode  array  before  and  after  lesion  formation  at 
normal  exposure  level  (2150  W/cm2  for  3  seconds).  The 
spectrograms  show  that  the  echoes  before  lesion  formation 
are  centered  at  I  MHz  with  small  2nd  harmonic  compo¬ 
nent  in  the  tissue  region  (at  100  mm).  The  spectrogram  of 
the  RF  data  from  the  same  direction  after  lesion  formation 
shows  a  strong  2nd  harmonic  component  starting  at  82  mm. 
This  is  consistent  with  the  shape  of  the  lesion  determined  by 
histological  evaluation  fl]  (tear-drop cross  section  with  the 
tip  near  the  geometric  center  and  the  base  near  the  front  of 
the  sample).  It  is  interesting  to  note  that,  in  addition  to  the 


increased  2nd  harmonic  generation,  the  echoes  from  the  le¬ 
sion  appear  to  have  wider  bandwidth  compared  to  echoes 
from  the  same  location  before  lesion  formation.  This  is 
reminiscent  of  the  behavior  of  ultrasound  contrast  agents 
(UCAs).  It  provides  indirect  evidence  that  bubble  oscilla¬ 
tion  occurs  during  and  after  lesion  formation  and  lingers  for 
tens  of  seconds  (in  vitro  depending  on  lesion  size). 

In  [1]  we  have  shown  images  of  single  shot  HIFU  le¬ 
sions  at  the  fundamental  and  SH  frequencies  of  the  dual¬ 
mode  transducer.  Both  imaging  modes  produced  satisfac¬ 
tory  images  of  the  lesion  in  these  cases,  with  the  SH  im¬ 
ages  providing  higher  contrast  of  lesion  echogenicity.  How¬ 
ever,  both  of  these  methods  suffer  when  imaging  extended 
lesions.  One  such  example  is  shown  in  Figure  4  before  (left) 
and  after  lesion  formation  at  the  fundamental  (top  pair),  SH 
(middle  pair)  and  using  the  quadratic  component  from  the 
SVF  (bottom  pair).  One  can  easily  observe  the  increase  in 
echogenicity  at  the  lesion  location  (90  -  100  mm  along  the 
axis  of  the  array)  for  all  three  modes.  However,  one  can  see 
that  the  dynamic  range  of  the  fundamental  image  is  limited 
by  the  speckle  pattern  from  tissue  surrounding  the  lesion 
while  the  SH  image  is  limited  by  the  beamforming  artifacts. 
The  lesion  echogenicity  in  this  case  is  nearly  17  dB  above 
the  normal  tissue  echogenicity  in  standard  echographic  im¬ 
age  (typical  contrast  is  more  like  5  -  7  dB).  The  SH  images 
improve  the  lesion  contrast  to  approximately  25  dB  com¬ 
pared  to  normal  tissue  (22  -  25  dB  typical),  but  beamform¬ 
ing  artifacts  compromise  the  definition  of  the  lesion  bound¬ 
aries  in  the  lateral  direction. 

Using  a  beamformed  echo  data  segment  along  the  main 
axis  of  the  array,  the  coefficients  of  linear  and  quadratic 
components  of  the  SVF  were  estimated  as  described  in  Sec¬ 
tion  3.2.  The  quadratic  components  of  the  beamformed  RF 
data  were  filtered  out  at  all  pixel  locations.  The  images  be¬ 
fore  and  after  lesion  formation  are  shown  in  Figure  4  at  40 
dB  dynamic  range.  One  can  see  quite  clearly  that  the  grating 
lobe  components  visible  in  the  2nd  harmonic  images  arc  ef¬ 
fectively  eliminated  in  the  quadratic  component  images  ob¬ 
tained  through  the  SVF.  In  addition,  the  speckle  components 
conspicuous  in  the  standard  echo  images  is  greatly  reduced. 
This  level  of  enhancement  is  typical  and  has  been  observed 
consistently  in  over  100  experiments  similar  to  the  one  de¬ 
scribed  in  this  paper.  The  reader  can  appreciate  that  the  le¬ 
sion  boundaries  are  well  defined  in  both  the  axial  and  lateral 
directions.  With  the  significant  increase  in  dynamic  range, 
one  can  see  that  both  detection  and  mapping  of  thermal  le¬ 
sions  is  significantly  facilitated  by  the  use  of  the  quadratic 
filter  based  on  the  SVF  model. 

5.  CONCLUSIONS 

Experimental  results  from  ex  vivo  tissue  samples  provide 
the  strongest  evidence  yet  that  thermal  lesions  exhibit  non¬ 
linear  behavior  as  a  propagation  medium.  Using  a  the  SVF 
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Fig.  2.  R F  data  and  spectrograms  before(left)  and  after 
(right)  normal  exposure. 
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Fig.  4.  Images  before  (left)  and  after  (right)  formation  of 
lesion.  Top:  standard  echo.  Middle:  second  harmonic.  Bot¬ 
tom:  Quadratic. 


Fig,  3.  RF  data  and  spectrograms  before  (left)  and  after 
(right)  over  exposure. 


model,  we  have  separated  the  linear  and  quadratic  compo¬ 
nents  in  the  beamformed  RF  echo  data  and  formed  images 
from  the  quadratic  components.  The  quadratic  component 
images  show  significant  enhancement  in  lesion  visualiza¬ 
tion  due  to: 

1  They  directly  exploit  the  nonlinear  nature  of  freshly 
formed  thermal  lesion  (possibly  due  to  formation  of 
microbubbles). 

2.  Quadratic  component  combines  both  low  frequency 
(close  to  dc)  and  harmonic  frequency  in  forming  non¬ 
linear  echoes.  This  simultaneously  reduces  speckle 
and  beamforming  artifacts  without  loss  in  spatial  res¬ 
olution. 

3.  The  quadratic  kernel  of  the  SVF  rejects  the  additive 
white  Gaussian  noise  components  which  significantly 
improves  the  SNR  of  the  imaging  system  and  enhances 
the  visualization  of  low  echogenicity  regions  in  the 
image. 
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ABSTRACT 

Nonlinear  ultrasonic  imaging  methods  (like  pulse  inversion1  and  quadratic  imaging 
based  second-order  Volterra  filter8)  are  used  in  visualization  of  lesion  formation 
in  freshly  excised  tissue.  Both  of  these  methods  are  more  sensitive  to  nonlinear 
echoes  (e.g.,  due  to  micro- bubbles)  than  standard  B-mode  imaging.  While  all  three 
methods  typically  show  increased  echogenicity  at.  the  lesion  location,  the  nonlinear 
methods  exhibit  more  localized  echo  enhancement  than  B-mode  imaging.  There¬ 
fore,  nonlinear  methods  are  potentially  better  suited  to  lesion  mapping  for  purposes 
of  image  guidance.  Quadratic  images  have  the  added  advantage  of  a  significant  in¬ 
crease  in  image  dynamic  range  and  noise  reduction  (a  major  limitation  of  pulse 
inversion  imaging).  The  results  shown  in  this  report  continue  to  support  the  hy¬ 
pothesis  that  micro-bubbles  play  an  important  role  of  lesion  formation.  In  this 
paper,  we  present  imaging  results  before  and  after  volumetric  lesion  formation  in 
ax  vivo  tissue.  The  results  illustrate  the  advantage  of  nonlinear  imaging  methods 
compared  to  conventional  B-scan  imaging  in  terms  of  accurate  mapping  of  lesion 
size  and  location. 

Keywords:  Noninvasive  surgery,  image  guidance,  treatment  monitoring,  ultrasonic 
imaging,  phased  arrays,  autoregressive  modelling. 


1.  INTRODUCTION 

We  have  used  a  dual-mode  array  described  in2  to  form  HIFU-induccd  thermal  le¬ 
sions  in  freshly  excised  degassed  tissue  under  a  variety  of  normal  exposure  and  over 
exposure  conditions.  Single-transmit  focus  images  were  collected  for  over  100  lesions 
before  and  after  lesion  formation.  These  images  have  consistently  shown  5  -  7  dB  en¬ 
hancement  in  the  echogenicity  from  the  lesion  location  in  the  standard  echographic 
images.  These  results  were  much  more  consistent  than  the  reported  “Hashes”  on  the 
B-scan  images  when  diagnostic  ultrasound  systems  are  used  to  monitor  HIFU  lesion 
formation.  Motivated  by  the  excellent  investigation  by  P.  P.  Lele  reported  in,5  we 
hypothesized  that  this  change  in  echogenicity  is  due  to  stable  inicrobubbles  that 
can  occur  even  at  low  insonation  levels.  Lele  found  that  subharmonic  emission  due 
to  microbubbles  showed  a  monotonic  increase  with  intensity  from  1.50  m\V/<;nr  to 
1500  Wr/crn2  without  a  distinct  threshold  for  emission  (measurements  done  in  vitro 
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and  in  two  at  2.7  and  1.8  MHz).  The  consistency  of  the  increase  in  echogenicity 
at.  t  le  lesion  may  be  explained  by  the  fact  that  the  microbubbles  may  already  be 
resonant  at  the  imaging  frequency  (same  as  the  therapeutic  HIFU  beam  when  the 
dual-mode  array  is  is  used),  perhaps  a  result  of  rectified  diffusion. 

The  standard  echographic  images  at  the  fundamental,  however,  offer  limited 
contrast  enhancement  due  to  the  speckle  phenomenon.  Therefore,  they  could  not 
provide  a  reliable  method  for  mapping  the  boundaries  of  HIFU-induced  lesions.  This 
lead  us  to  try  to  exploit  the  nonlinear  nature  of  the  microbubbles  to  enhance  the 
visualization  and  mapping  of  thermal  lesions.  The  idea  was  that,  if  microbubbles  are 
indeed  present  at  the  lesion  location,  they  will  generate  nonlinear  echoes  that  may 
be  better  suited  for  mapping.  We  have  initially  investigated  second  harmonic  (SH) 
imaging  as  a  means  of  enhancing  the  lesion  contrast  for  improving  the  visualization 
of  these  images.  SH  images  of  thermal  lesions  have  shown  increase  in  the  contrast  on 
the  order  of  22  -  25  dB.  but  with  decreased  dynamic  range  of  the  resulting  images.2 
A  post-beamforining  nonlinear  compounding  algorithm  was  shown  to  improve  the 
contrast  in  lesion  echogenicity  to  30  -  35  dB  without  loss  in  dynamic  range.7  This 
was  achieved  by  compounding  the  fundamental  and  the  SH  images  using  spatial 
compounding  functions  based  on  the  receive  beamforming  characteristics  of  the 
dual-mode  array  at  the  fundamental  and  the  SH  frequencies.  In  this  paper,  we  use 
a  commercial  imaging  scanner  with  modifications  to  allow  pulse  inversion  imaging  in 
addition  to  standard  B-mode  imaging  for  the  visualization  of  freshly  excised  tissue 
before,  during,  and  after  the  formation  of  HIFU  lesions. 

2.  NONLINEAR  IMAGING  METHODS 

2.1.  Pulse  Inversion  Imaging 

This  method  was  recently  introduced  by  Burns  and  coworkers6  for  enhancing  con¬ 
trast  echoes  in  contrast-assisted  pulse-echo  imaging.  The  basic  implementation  of 
this  method  entails  the  use  of  two  pulses  per  image  line.  These  pulses  are  carefully 
designed  such  that.  p2{t)  ~  -p\[t  -  TL),  where  T/,  is  some  appropriate  delay  (on 
the  order  of  the  pulse-echo  time  from  the  maximum  depth  of  interest).  Summing 
the  edioes  resulting  from  the  two  pulses  eliminates  the  odd-harmonic  components 
from  the  echo  signal  (including  the  fundamental)  while  doubling  the  even-harmonic 
(mostly  second)  components.  This  method  currently  represents  the  leading  ap¬ 
proach  for  contrast-agent  imaging,  especially  with  low  concentration  and/or  very 
low'  transmit  signals  to  minimize  the  generation  of  tissue  nonlinearity. 

2.2.  Quadratic  Imaging  Based  On  SVF 

We  have  recently  developed  a  now'  nonlinear  imaging  system  based  on  the  SVF.s 
This  lias  a  number  of  advantages  when  compared  to  PI  imaging  (e.g.  requires  only 
a  single  pulse  per  line  and  increased  dynamic  range). 
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2.2.1.  Second-Order  Volterra  Model 

Results  from  have  shown  the  validity  of  a  second-order  Volterra  filer  as  a  model 
for  pulse-echo  ultrasound  imaging  data  from  tissue  mimicking  media.  In  this  sec¬ 
tion,  the  decomposition  of  received  echo,  i.e.,  output  sequences  only,  into  linear 
and  quadratic  components  by  using  least-squares  approach  of  second-order  Volterra 
model  will  be  considered  and  the  detail  of  algorithm  implementation  to  pulse-echo 
ultrasound  imaging  will  be  stated. 

2.2.2.  Signal  Separation  Model 

The  algorithm  described  in  this  section  is  adapted  from.3  The  response  of  a 
quadraticallv  nonlinear  system,  y{n-r  1).  can  be  predicted  by  a  second- order  Volterra 
model  of  past  m.  values  as  follows: 

=  yL(n  +  1)  +  yQ(n  +  1) 

77i-  J  rfi —  1 

=  y(n  ~  OM*)  +  ^  i/(n  -j)y(n  -  k)hQ(j%k)  +  e(»).  (1) 

1=0 .  j= 0  k-j 

wheie  hi{i)  is  linear  filter  coefficients,  hq(j,  k)  represents  quadratic  filter  coefficients 
and  e(n)  is  a  modelling  error  and/or  a  measurement  noise  which  is  assumed  to  be  an 
independent,  identically  distributed(i.i.d)  random  variable  with  zero  mean.  That  is, 
if  the  model  coefficients  are  known,  the  echo  signal  can  be  decomposed  into  linear 
and  quadratic  components.  The  latter  can  be  expected  to  better  represent  the 
quadratic  response  of  the  system  than,  say.  the  second  harmonic  component.  The 
model  coefficients  can  be  obtained  by  setting  the  linear  and  quadratic  prediction 
problem  in  Equation  1  to  form  a  set  of  linear  equations.  Recognizing  that  the 
output  is  linear  in  terms  of  the  (unknown)  model  coefficients,  one  obtains  a  matrix 
equation  of  the  form: 

f  =  Gh  +  f,  (2) 

whore  the  vector  f,  the  matrix  G  and  the  error  vector  t  are 

f  =  [if(n  +  l),jK»  +  2) . y{n  +  L)]T 

G  =  [y(n).y(n+  l) . y(n  +  L-l)]T 

€  =  [e(n),5(n  +  l) . s(n  +■  L  -  1)]T. 

where  the  data  vector,  y,  is  given  by: 

y(n)  =  [2/(n),  y[n  -  1),  y(n  -  2), . . .  ,y(n  -  m  +  1). 
y"{n).y{ri)y(n  -  1) . i f{n  -  m  +  l)]r 

ami  the  filter  coefficient  vector,  h.  is  given  by: 

h  =  MbMl),/iL<2) . hL(rn-  1), 

0),  Iiq( 0, 1) . hQ(m  -  l,m  -  l)]7  . 

Hie'  details  of  the  solution  for  the  coefficients  of  the  SVF  model  can  be  found  in.4 
Briefly,  a  minimum-norm  least-squares  solution  of  (2)  is  obtained  using  truncated 
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siiuiular  value  decomposition.  TSVD.  To  assess  the  performance  of  the  signal  separa¬ 
tion  mode!  in  enhancing  the  lesion  visualization,  we  compute  the  contrast- to-t issue 
ratio: 

CTB=l0be'»(Ki) 

where  Wyqch  and  Ijy^rlb  are  the  l2  norms  of  the  quadratic  components  from  the 
lesion  and  normal  tissue  regions,  respectively.  These  regions  are  easily  identified 
under  various  imaging  conditions.  For  instance,  for  the  application  described  in 
this  paper,  the  contrast  region  is  the  expected  location  of  the  thermal  lesion  (often 
visible  on  the  standard  echographie  image). 

2.2.3,  Implementation 

The  coefficients  of  the  linear  and  quadratic  components  of  the  SVF  model  are 
obtained  from  the  beamformed  RF  data.  The  implementation  steps  are  as  follows: 

1.  Using  the  standard  echographie  image,  select  a  beamformed  RF  data  segment 
from  the  expected  lesion  location. 

2.  Form  the  linear  systems  of  equations  according  to  (2). 

3.  Define  contrast  region  (within  the  lesion)  and  normal  tissue  region  for  the 
computation  of  the  mean-square  error  MSE  and  CTR. 

4.  Solve  systems  of  linear  equations  by  using  TSVD  regularization  method. 

o.  Apply  second-order  Volterra  filter  to  the  beamformed  RF  data  throughout 
the  pulse-echo  ultrasound  image. 

6.  The  quadratic  component  from  the  SVF  can  be  displayed  as  a  separate  image 
or  appropriately  compounded  with  the  linear  component. 

3.  EXPERIMENTAL  SETUP 

Figure  1  shows  a  simple  arrangement  for  the  formation  of  HIFU  lesions  in  freshly 
excised  and  degassed  portin''  livers  samples.  The  therapy  transducer  is  a  1.3  MHz 
single-element  spherical-shell  transducer  with  a  radius  of  curvature  equal  to  its  di¬ 
ameter  and  equal  to  63.5  mm  (Etalon.  Lizton.  Indiana).  The  transducer  is  fixed  to 
the  hack  of  a  small  tank  as  shown  and  driven  by  a  power  amplifier  ( ENT.  Rochester. 

)  and  a  programmable  function  generator.  This  assembly  can  be  used  in  gen¬ 
erating  a  variety  of  amplitude-modulated  HIFU  bursts  from  tens  of  milliseconds  to 
several  seconds  long  and  intensities  up  to  3000  W/cm2  (conservative  estimate). 


Therapy  Transducer  (Thx  Tx) 


Imaging  Transducer  (lx  Tx) 
Acoustic  Window 


Figure  1.  Experiment  setup  used  for  formation  of  HIFU  lesions  using  a  single- 
element,  transducer  attached  to  the  of  the  tank  with  an  imaging  probe  monitoring 
a  cross  section  of  the  lesion  through  an  acoustic  window. 


3.1.  Real-time  Imaging 

Real-time  imaging  is  performed  using  a  modified  Technos  MP  system  from  ESAOTE. 
Genoa.  Italy.  The  system  is  modified  to  allow  imaging  in  puise  inversion  mode  in 
addition  to  normal  B-mode  imaging.  In  addition,  a  hardware  module  for  capturing 
high-quality  beamformed  RF  data  allows  us  to  capture  and  upload  up  to  60  seconds 
of  lull  frame  data  with  a  specified  frame  rate.  A  CA  421  convex  probe  was  used  in 
acquiring  image  data  for  this  paper.  Image  data  was  acquired  in  pulse  inversion  (PI) 
mode  with  a  2-cyele  transmit  pulse  centered  at  1.57  MHz.  The  imaging  transducer 
was  aligned  so  that  the  image  plane  is  x  -  y  plane  (to  allow  imaging  a  cross  section 
of  the  lesion).  However,  imaging  in  the  x  -  z  plane  was  often  used,  especially  for 
imaging  discrete  (single-shot)  lesions  (to  allow  imaging  the  lesion  along  the  axis  of 
the  therapeutic  array). 

3.2.  Lesion  Formation 

Volumetric  lesions  were  formed  by  continuously  driving  the  therapy  transducer  while 
moving  the  sample  holder  in  a  raster  scan  with  the  stepper  motor.  The  motor  speed, 
input  voltage  profile,  and  line  sparing  in  the  raster  were  changed  to  produce  a  variety 
of  volumetric  lesions  that  can  be  characterized  as  over  exposure ,  normal  exposure 
or  under  erasure. 
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4.  RESULTS  AND  DISCUSSION 

The  experimental  setup  shown  in  Figure  1  was  used  in  obtaining  images  of  ex  vivo 
tissue  samples  before  and  after  thermal  lesion  formation  with  the  spherical  therapy 
transducer.  This  setup  was  designed  to  allow  comparisons  with  our  a  dual-mode 
64-element  ultrasound  phased  array  operating  at  1  MHz  described  in.2  Results 
shown  in-  confirmed  that  echoes  from  the  lesion  location  exhibited  increased  levels 
of  second  harmonic  generation  even  at  normal  exposure  conditions.  In  this  paper, 
we  show  imaging  results  from  continuous  exposure  at  1800  W/cm2  while  the  tissue 
is  moved  in  a  raster  scan  10  x  10  mm  with  raster  lines  parallel  to  the  face  of  the 
imaging  transducer  (1  mm  spacing).  The  motor  speed  was  set  at  2  nnn/s  tints 
giving  an  exposure  time  of  roughly  1  second  at  one  location  in  the  focal  plane  of 
the  therapeutic  transducer,  but  larger  in  the  prefocal  region.  Figure  2  shows  a  cross 
section  of  the  tissue  after  lesion  formation  with  conspicuous  thermal  lesion  in  the 
upper  left  corner.  The  cross  section  conforms  approximately  to  the  programmed 
pattern  with  distortions  most  probably  due  to  deformation  of  tissue  as  it  was  cut 
without  fixing.  The  tissue  is  positioned  so  that  the  therapeutic  transducer  is  aimed 
upwards  into  the  imaging  plane  and  the  imaging  transducer  is  imaging  downwards. 
Gross  measurements  of  the  cross  section  of  the  lesion  shows  a  lateral  extent  of  10 
rum  and  an  axial  extent  of  12  (with  respect  to  the  imaging  transducer). 

Figure  3  shows  B-mode  (left),  PI  (center),  and  quadratic  (right)  images  before 
(top)  and  after  lesion  formation.  All  imaging  modes  show  enhanced  echogenicity  at 
the  lesion  location  at  an  axial  distance  of  145  mm  from  the  imaging  transducer.  The 
PI  and  quadratic  images  show  a  smaller  size  hyperechoic  region  with  lateral  extent 
of  about  12  mm  compared  to  13  nun  in  the  B-mode  image.  The  axial  extent  of  the 
hyperechoic  region  is  11  mm  in  the  PI  image  and  10  mm  in  the  quadratic  image. 
It.  was  difficult  to  determine  this  value  in  the  B-mode  image  due  to  the  interference 
of  the  speckle  pattern  in  tissue.  This  result  is  typical,  i.e.  the  PI  data  gives  a  more 
accurate  map  of  the  lesion  than  B-mode  data. 

Comparing  the  quadratic  and  B-mode  images,  one  can  see  the  speckle  compo¬ 
nents  conspicuous  in  the  standard  echo  images  is  greatly  reduced.  This  level  of 
enhancement  is  typical  and  has  been  observed  consistently  in  over  100  experiments 
similar  to  the  one  described  in.2  The  reader  can  appreciate  that  the  lesion  bound¬ 
aries  are  well  defined  in  both  the  axial  and  lateral  directions.  With  the  significant 
increase  in  dynamic  range,  one  can  see  that  both  detection  and  mapping  of  thermal 
lesions  is  significantly  facilitated  by  the  use  of  the  quadratic  filter  based  on  the  SVF 
model.  It  is  also  interesting  to  not  the  remarkable  similarity  between  the  PI  and 
quadratic  images,  except  that  the  latter  has  more  than  twice  the  dynamic  range. 
These  results  are  also  typical  and  represent  the  potential  advantage  of  the  quadratic 
processing  versus  pulse  inversion.  The  result,  can  be  understood  by  keeping  in  mind 
that  both  methods  are  sensitive  to  the  harmonic  content  of  the  RF  data,  but  they 
differ  in  their  respective  degrees  of  suppressing  the  noise  component;  the  quadratic 
signal  being  far  superior  in  this  respect. 


Figure  2.  A  section  of  a  volumetric  lesion  in  an  ex  vivo  porcine  liver  sample. 


-  ,  5.  CONCLUSIONS 

experimental  results  from  ex  vivo  tissue  samples  provide  compelling  evidence  that 
thermal  lesions  exhibit  nonlinear  behavior  as  a  propagation  medium.  Nonlinear 
imaging  methods  were  used  to  separate  the  linear  and  nonlinear  components  in 
the  beamlormed  RF  echo  data.  PI  images  confirmed  the  presence  of  strong  2nd 
harmonic  component  in  echoes  from  lesion  location.  The  quadratic  component 

images  (obtained  through  SVF)  show  significant  enhancement  in  lesion  visualization 
due  to: 

1.  They  directly  exploit  the  nonlinear  nature  of  freshly  formed  thermal  lesion 
(possibly  due  to  formation  of  microbubbles). 

2.  Quadratic  component  combines  both  low  frequency  (close  to  dc)  and  harmonic 
frequency  in  forming  nonlinear  echoes.  This  simultaneously  reduces  speckle 
and  beamforming  artifacts  without  loss  in  spatial  resolution. 

3  Thc  Quadratic,  kernel  of  the  SVF  rejects  the  additive  white  Gaussian  noise 
components  which  significantly  improves  the  SNR  of  the  imaging  system  and 
enhances  the  visualization  of  low  echogenicity  regions  in  the  image. 

Efforts  to  further  characterize  echo  signals  from  the  different  imaging  methods  and 
correlate  them  to  the  state  of  the  tissue  are  currently  underway. 
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Figure  3.  B-mode  (left),  pulse  inversion  (center),  and  Quadratic  images  before 
(upper)  and  immediately  after  (lower)  formation  of  volumetric  lesion.  Dynamic 
range:  B-mode  60dB,  PI  30  <1B.  Quadratic  75  dB. 
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Abstract — We  present  a  new  algorithm  for  deriving  a 
second-order  Volterra  filter  (S  VF)  capable  of  separating  lin¬ 
ear  and  quadratic  components  from  echo  signals.  Images 
based  on  the  quadratic  components  are  shown  to  provide 
contrast  enhancement  between  tissue  and  ultrasound  con¬ 
trast  agents  (UCAs)  without  loss  in  spatial  resolution.  It 
is  also  shown  that  the  quadratic  images  preserve  the  low 
scattering  regions  due  to  their  high  dynamic  range  when 
compared  with  standard  B-mode  or  harmonic  images.  A 
robust  algorithm  for  deriving  the  filter  has  been  developed 
and  tested  on  real-time  imaging  data  from  contrast  and 
tissue-mimicking  media.  Illustrative  examples  from  image 
targets  containing  contrast  agent  and  tissue-mimicking  me¬ 
dia  are  presented  and  discussed.  Quantitative  assessment 
of  the  contrast  enhancement  is  performed  on  both  the  RF 
data  and  the  envelope-detected  log-compressed  image  data. 
It  is  shown  that  the  quadratic  images  offer  levels  of  enhance¬ 
ment  comparable  or  exceeding  those  from  harmonic  filters 
while  maintaining  the  visibility  of  low  scattering  regions  of 
the  image. 


1.  Introduction 

INc  :RUASlNg  interest  in  extending  the  capabilities  of  ul¬ 
trasound  imaging  by  utilizing  ultrasound  contrast  agents 
(UCAs)  has  heightened  the  need  for  more  suitable  imaging 
techniques.  In  standard  B-mode  imaging.  UCAs  increase 
the  echogenicity  from  perfused  tissues  [1 1.  This  results  in 
improved  endocaidial  border  detection  in  left  ventricular 
opacification,  which  leads  to  a  better  analysis  of  wall  mo¬ 
tion  abnormalities  [2].  Nevertheless,  in  the  myocardium 
where  the  ratio  of  blood  volume  to  tissue  is  quite  low 
( approximately  10%  [3]),  the  backscatter  from  the  small 
number  of  microbubbles  in  vessels  can  be  dominated  by 
echoes  from  surrounding  tissue.  In  this  occurrence,  stan¬ 
dard  B-mode  imaging  offers  inferior  UCA  detectability  in 
the  presence  of  tissue,  stated  as  agent-to- tissue  ratio  [4]. 
hi  order  to  increase  the  sensitivity  of  UCA  detections,  var¬ 
ious  new  imaging  techniques  [5],  [6]  have  been  developed 
by  employing  some  specific  acoustic  signatures  of  UCAs. 
such  as  nonlinear  and  transient  scattering. 

Imaging  techniques  based  on  nonlinear  oscillations  have  > 
been  designed  for  separating  and  enhancing  nonlinear 
UCA  echoes  from  a  specified  region  of  interest  within  the 
imaging  field,  including  second  harmonic  (SH)  B-mode 
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imaging  and  pulse  inversion  (PI)  Doppler  imaging  :7 
The  SH  imaging  employs  a  fundamental  frequency  tiwis 
init  pulse  and  produces  images  from  the  second  harmnnk 
component  of  received  echoes  bv  using  a  second  hannonn 
bandpass  filter  (BPF)  to  remove  the  fundamental  fre¬ 
quency.  In  order  to  increase  UCA  detection  seusitivin  in 
the  limited  transducer  bandwidth  condition,  spectral  over¬ 
lap  between  fundamental  and  second  harmonic  part,*  need 
to  be  minimized  by  transmitting  narrow- hand  pulses  re¬ 
sulting  in  an  inherent  tradeoff  between  contras i  and  spa¬ 
tial  resolution. 

In  PI  imaging,  a  sequence  of  two  inverted  acoustic 
pulses  with  appropriate  delay  is  transmitted  into  tissue. 
Images  are  produced  by  summing  the  corresponding  two 
backscattered  signals.  In  the  absence  of  tissue  motion,  ihe 
resulting  sum  can  be  shown  to  contain  only  even  harmonies 
of  the  nonlinear  echoes  [7].  The  PI  imaging  overcomes  t  lie 
tradeoff  between  contrast  and  spatial  resolution  because  it 
utilizes  the  entire  bandwidth  of  the  backs<  altered  signals 
[7].  As  a  result,  superior  spatial  resolution  can  be  achieved 
when  compared  with  SH  imaging.  Moreover,  it  has  been 
shown  that  PI  imaging  can  be  operated  in  a  continuous 
imaging  inode  with  low  mechanical  indices  (Mis)  [8!  The 
PI  imaging  is  sensitive  to  tissue  motion  because  it  is  a 
multiple  pulse  technique:  therefore,  PI  detection  is  com¬ 
bined  with  Doppler  detection  leading  to  a  new  technique 
called  PI  Doppler.  The  PI  Doppler  utilizes  the  advantages 
from  both  detection  schemes  and  circumvents  the  tissue 
motion  problem  [7].  Nevertheless,  an  inherent  mult  ipttlse 
technique  of  PI  imaging  results  in  the  reduction  of  imaging 
frame  rates. 

In  order  to  detect  backscattered  signals  due  to  UCAs  in 
pulse-echo  ultrasound  imaging  with  single  transmit  pulse, 
per  line  and  overcome  the  tradeoff  between  contrast  and 
spatial  resolution,  we  have  developed  a  new  imaging  tech¬ 
nique  based  on  the  Volterra  filter  [9]  The  Volterra  filler 
is  a  dynamic  filter  that  operates  in  parallel  on  the  lin¬ 
ear.  quadratic,  cubic,  etc.,  signal  components  to  produce 
its  output.  It  has  a  discrete  convolutional  form  with  finite 
memory  that*  is  commonly  used  in  nonlinear  DSP  [AU: 
Define  DSP  on  first  use.]  applications  [10].  While  the 
Volterra  filter  output  is^nonlinear  in  terms  of  the  input 
data,  it  is  linear  in  terms  of  its  coefficients.  As  a  result, 
linear  processing  can  be  applied  to  identify'  Volterra  ker¬ 
nels  (i.e..  filter  coefficients).  The  identification  of  Yohetra 
kernels  has  been  demonstrated  in  several  application*  ‘Hi:. 
[31].  For  example,  the  linear  and  quadratic  Volieria  ker¬ 
nels  identified  by  an  adaptive  filtering  aigoiiriini  i >«-;*■?!  uii 
recursive  least-squares  approach  of  a  s»erond-oi  u-t  Y«  coyra 
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model  are  applied  to  separate  the  linear  and  quadratic  re* 
spouses  of  a  tension  leg  platform  [12]. 

In  this  paper,  we  present  post  beam  forming  nonlinear 
filter  based  on  the  second-order  Volterra  filter  (SVF).  The 
filter  is  capable  of  separating  linear  and  quadratic  compo¬ 
nents  from  UGA-backscattered  signals.  Filter  coefficients 
can  be  identified  by  forming  a  system  of  linear  equations 
from  a  beamformed  RF  data  segment  in  the  UCA  or  tis¬ 
sue  region  and  obtained  by  solving  a  minimum-norm  least- 
squares  (MNLS)  problem.  It  is  shown  that  the  system  iden¬ 
tification  approach  to  a  nonlinear  ultrasonic  imaging  is  ro¬ 
bust.  and  can  be  applied  to  form  nonlinear  images  through¬ 
out  the  imaging  field,  i.e..  not  confined  to  the  region  where 
t  he  filter  coefficients  were  derived.  Images  produced  from 
quadratic  output  of  the  SVF  model  exhibit  excellent  con¬ 
trast.  enhancement  without  loss  in  spatial  resolution.  In 
addition,  quadratic  images  have  increased  dynamic  range 
compared  to  standard  B-niode  and  harmonic  images  allow¬ 
ing  the  preservation  of  image  features  along  with  contrast 
improvement. 


II.  Theory 

For  simplicity  and  without  loss  of  generality’,  we  illus¬ 
trate  the  nonlinear  postbeam  forming  algorithm  using  the 
SVF  model.  The  algorithm  described  in  this  section  for  de¬ 
riving  the  filter  coefficients  of  the  SVF  extends  to  higher 
oider  in  a  straightforward  manner.  Initial  experience  with 
ihis  model  indicate  that  it  is  largely  sufficient  for  tissue  re¬ 
sponse.  Furthermore,  the  SVF  is  computationally  feasible 
for  real-time  application  with  today's  technology,  which 
makes  it.  more  attractive  from  the  implementation  point 
of  view. 

A.  Second- Order  Volterra  Model 

RcsuIts  from  [9]  have  shown  the  validity  of  a  SVF  as 
a  model  for  pulse-echo  ultrasound  pulse-echo  data  from 
tissue  mimicking  media.  An  input-output  system  identifi- 
cai  ion  approach  was  used  to  estimate  the  coefficients  of 
the  linear  and  quadratic  components  of  the  SVF  in  the 
frequency  domain.  However,  while  the  system  identifica¬ 
tion  st  udy  was  necessary  to  establish  the  applicability  of 
SVF  to  ultrasound  pulse-echo  data,  it  is  not  useful  for 
imaging  purposes  as  it  requires  access  to  both  the  input 
and  the  echo  data  from  distinct*  scat-1  erers  in  the  tissue- 
mimicking  media.  An  appropriate  approach  for  imaging 
operates  on  the  beamformed  RF  data  to  separate  the  lin¬ 
ear  and- quadratic  components  regardless  of  the  input  .  This 
signal  separation  approach  allows  us  to  extract  the  linear 
and  quadratic  signal  components  from  the  beamformed 
dat  a  to  form  linear  and/or  quadratic  images  separately  or 
compounded.  In  this  paper,  we  emphasize  the  quadratic 
image*  obtained  from  the  SVF  model  based  on  a  linear 
and  quadratic  prediction  model  of  boa  informer  output  de¬ 
scribed  in  the  remainder  of  this  section. 


u{n) 

Beamformer  - 


*  O' ,<  /. k  )  ►  ?!  •a) 

Fig.  1.  Separation  of  beamformed  RF  data  into  linear  and  quadiatir 
components  using  the  SVF. 

L  Signal  Separation  Model:  Fig.  1  shows  a  simple  block 
diagram  of  the  imaging  system  based  on  SVF.  The  SVF 
operates  on  the  beamformer  output  to  product*  the  linear 
and  quadratic  components.  uiJn)  and  uq{v),  respectively. 
Estimates  of  the  total  beamformer  output  can  be  obtained 
from  these  components  simply  by  adding  them 

u(n)  =  uL(n)  4  uo[ri ),  0  ) 

where  «(«),  ui(n ),  and  uq(v)  are  the  total,  linear,  and 
quadratic  estimations,  respectively.  The1  separation  of  the 
linear  and  quadratic  components  can  be  achieved  once 
t  he  coefficients  of  the  kernels,  hL(i)  (linear ).  and  ) 

(quadratic)  are  found.  In  the  following  subsection,  we  de¬ 
scribe  a  MNLS  approach  for  determining  these  coefficients. 

2.  The  MNLS  Estimation  of  SVF  Coefficients:  The  re¬ 
sponse  of  a  quadratically  non  linear  system  with  memory. 
din  4-  1)  can  be  predicted  by  a  (discrete)  second-order 
Volterra  model  operating  on  the  in  past  samples  as  fol¬ 
lows 

m  —  ] 

u.(ti  +  1)  =  w(n  -  i)hi,li) 
m  —  l  m  —  1 

4*  u^n  “  jMri  ~  k)hQ(j.k).  (2) 

l;=j 

where  hL(i)  is  linear  filter  coefficients,  and  h.Q(j.h)  repre¬ 
sents  quadratic  filter  coefficients.  Note  that  while  u(n  4  1 ) 
is  nonlinear  with  respect  to  the  beamformed  data,  it  is 
linear  with  respect,  to  the  coefficients  of  the  linear  and 
quadratic  kernels  of  the  SVF.  Recognizing  this  fan.  one 
can  rewrite  (2)  in  vector  .form 

u(n  4- 1 )  =  ur(?Mh.  .  c3i 

where  the  data  vector.  u(r/),  is  defined  at  sample  n  as 

b- 

u(ti)  =  [u(n),  v{n  -  1),  u(n  -  2) . u{n  -  tn  4  1  i. 

u2(7?).  u(n)u(n  -  1) . u2hi  -  ni  -  1  )];  . 

and  the  filter  coefficient  vector,  h.  can  lx*  expressed  as 

h  =  7,,/.(CH.  hi[  1).  hi{2 ).. ....  hjfrn  -  1 !. 

llQiO.0).  h-Q  (0, 1 ). .  .  I;q  { v )  -  ] .  vt  —  1  i;7  : 
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where  m  is  t.he  system  order  and  superscript  r  denotes  the  B.  Regularization 
transpose.  The  total  number  of  independent  filter  coeffi¬ 
cients,  A',  is ’equal  to  [m2  +  3m)/2  assuming  a  symmetri-  The  SVD  of  G  forms  a  basis  for  regularization  bv  up- 

cal  quadratic  kernel  (i.e.,  hQ(j.k)  =  hQ(k.j)).  Similarly,  ,  propriate  selection  of  singular  modes  that  enhance  the  ro- 

u{n  +  -).  uln  -f  3), . . . ,  u(n  +  M)  can  be  represented  in  the  constructed  image  in  some  sense.  There  are  a  number  of 

orni  of  (3)  and  expressed  in  the  matrix  form  approaches  for  regularization  of  (8),  including  single  pa 

f  =  GIT  f4)  rameter  aaid  rank  reduction  regularization  ;14j  Th«-  lat¬ 

ter.  sometimes  referred  to  as  the  truncated  singular  value 
where  the  vector  f  and  the  matrix  G  are  defined  as  decomposition  (TSVD).  produces  a  solution  bv  trnncat 

f  =  +  u(n  j.  2),  u(ri  +  in§  the  number  of  singular  modes  of  G  with  the  smallest 

‘  .  ’  singular  values  below  a  certain  threshold.  The  Inn  order 

TSVD  solution  is  given  by 


G  -  :u(n.J,  u(n  +  1) . u(n+  M  -  l)]r. 

where  M  is  the  nimxber  of  linear  equations  (observations). 
The  linear  and  quadratic  filter  coefficients  can  be  esti¬ 
mated  by  seeking  an  appropriate  solution  of  (4).  Well- 
known  solutions  to  (4)  are  the  least  squares  (LS)  solu¬ 
tion  for  the  overdetermined  case  (more  constraints  than 
unknowns)  and  the  minimum  norm  (MN)  solution  in  the 
underdetermined  case  (less  constraints  than  unknowns)  A 
MNLS  solution  can  be  obtained  by 


Hmnls  =  G‘f.  (5) 


when*  G'  is  a  generalized  inverse  [13].  This  solution  ap¬ 
plies  to  both  the  overdetermined  and  underdetermined 
cases  as  u  accounts  for  the  effective  rank  of  I  he  matrix 
G.  In  general.  G  could  be  rank  deficient,  i.e..  has  rank 
r  S  n:in{.1/.  A  }.  When  this  is  the  case,  there  are  infinite 
number  of  solutions  to  (4)  that  produce  the  same  LS  er¬ 
ror,  |Jf  -  GhI5||T  The  MNLS  is  the  unique  solution  to 
(4)  with  minimum  norm,  i.e.,  ||hMNLS|i2  <  |]h^||2  VhL5. 
The  singula]  value  decomposition  (SVD)  of  G  is  given  by 

G  =  UEVr  ■ 


and 


=  D<*u,vf.  ■ 

m  1 


(6) 


•ti( 


where  the  truncation  parameter  k  <  r.  also  known  a>  the 
rank  of  the  approximation,  is  the  number  of  singular  n;«».  ies 
used  to  compute  the  estimate. 

The  regularization  is  guided  by  the  mean  square  error 
criterion,  approximated  by 


E{h)  -  10  log J0  I 


f-Gh,  |12  ' 
M 


i  10) 


where  ||  *  ||2  is  the  lo  norm. 

In  the  context  of  the  linear  and  quadratic  prediction  ap¬ 
proach  taken  in  this  paper,  the  MSE  decreases  nmnotom- 
cally  with  k.  A  criterion  for  choosing  an  appropriate  value 
of  k  is  needed.  In  the  context  of  contrast -agent  imaging, 
an  obvious  criterion  is  the  contrast- to-tissue  ratio  «CTK  i 


where_Fc  is  the  average  power  of  signals  in  a  UCA  region, 
and  Ft  is  the  average  power  of  signals  in  a  tissue  region 
The  average  power  of  signals  in  a  given  region.  P.  can  be 
expressed  as 


G’  =  VEtU7’ 

=E-vjUr. 


7=1 


(7) 


where  E  is  a  M  x  .V  diagonal  matrix  with  singular  val- 
ues  a,  >  a-  >  ...  >  c Jr  >  crT+}  =  . . .  =  <rr  =  0  (p  = 
:nin{;lL.Y}).  The  matrices  U  (A!  x  M)  and  V(Ar  x  N) 
are  formed  front  the  columns  {u,}££3  and  {v7};l3.  which 
are  the  orthogonal  eigenvectors  of  GGr  and  G7  G.  respec¬ 
tively  i!3j.  Using  (7).  the  MNLS  solution  to  (5)  is  then 
given  by 


h\lNLS  =  ^  V;  . 

'  rr 


(8) 


7=1 


Xou*  That,  if  G  is  full  rank,  the  MNLS  is  equivalent  to  the 
LS  solution  in  the  over  determined  case  and  to  the  MN  solu¬ 
tion  in  rhehniderdetonnmed  case.  For  example,  for  M  <  .V 
! underdetermined).  hush$  =  hji/.v  =  GT(GGr)“!f. 


,7=1  ./=] 
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where  x,j  is  the  signal  in  that  region.  Note  that  the  CTR 
used  in  this  algorithm  is  determined  from  quadratic  com¬ 
ponents  of  the  SVF  model.  The  CTR  provides  a  meaning¬ 
ful  stopping  criterion  for  TSVD  since  Elk)  is  rnouotoni- 
calh*  decreasing  or  nonincreasing  with  k 

One  may  use  a  variety  of  methods  to  obtain  a  regular¬ 
ized  solution  to  the  set  of  equations  involving  the.  Volirna 
kernels.  Examples  are  linear  or  nonlinear  programming  m 
cases  when  certain  constraints  on  the  filter  coefficients  arc- 
known  to  apply  (e.g..  positivity)  or  penalized  maximum 
likelihood  when  known  statistical  properties  of  ihe  da: a 
can  be  incorporated.  A  powerful  approach  seeks  a  solution 
to  a  constrained  optimization  problem  of  the  form 


min  Rr  subject  to  Gh  =  f. 

hQ 
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where  R2  is  an  appropriately  chosen  quadratic  ratio  to  be' 
minimized  (e.g..  inverse  of  CTR)  For  example.  R2  can  be 
chosen  to  reflect  the  inverse  of  the  CTR  before  log  com¬ 
pression.  In  general,  it  can  be  any  quantity  that  depends 
nri  (he  solution  hQ  that  is  to  be  minimized.  This  leads  to 
a  modifier!  form  of  the  TSVD  given  here 

=  '  (14) 

»=1  Ui  T  i*'i 


Srart 


From  a  Standard  Hchographic  Image  Define 
A  Comrast  Region 
A  Tissue  Region 

A  Segment  of  the  Beam  formed  kf  da:.*« 
System  Orders,  w  fioin  F\a  Q 

-  .  T. 


wheie  :  is  an  appropriately  chosen  threshold,  and  R;  is  HRni  a ’SwanoT 

quadratic  ratio  resulting  from  the  quadratic  kernel  ob-  Linear  Equations 

t.ajned  from  the  fth  singular  mode.  Regardless  of  the  reg-  . f=Cli 

ularizat.ion  procedure,  the  fundamental  result  here  is  the  . . .  V- . 

use  of  the  linear/quadratic  prediction  to  obtain  a  set  of  in-  Compute  svu  of  G 

dependent  equations  that  can  l>e  solved  robustly  to  obtain  -  .  - 

the  Volterra  filter  kernels.  r 


C  Quadratic  Images 

Quadratic  images  are  obtained  from  quadratic  compo¬ 
nents  of  the  second-order  Volterra  model.  The  coefficients 
of  the  SVF  are  derived  from  the  beamformed  RF  data 
taken  from  a  representative  region  on  a  standard  B-mode 
image  Details  of  the  algorithm  to  produce  the  quadratic 
image  are  as  follows. 

From  the  standard  B-niode  image,  a  UCA  region  arid 
a  tissue  region  are  defined  for  the  CTR  computation.  The 
definition  of  the  CTR  reference  regions  depends  on  the 
imaging  target  are  described  in  Section  II1-C.  In  general, 
we  try  to  find  regions  at  the  same  depth  and  with  the 
same  l>eam  angle  with  respect  to  the  axis  of  the  imaging 
array.  Furthermore,  whenever  possible,  we  chose  multiple 
overlapping  subregions  to  obtain  multiple  CTR  values  at 
different  depths. 

Once  the  CTR  reference  regions  are  defined,  a  segment 
of  RF  data  from  an  axial  line  is  selected  to  form  a  sys¬ 
tem  of  linear  equations  according  to  (4).  This  segment 
can  be  selected  from  the  tissue  or  the  UCA  region  as  long 

the  appropriate  regularization  of  the  MNLS  solution  is 
sou glit.  For  example,  when  TSVD  is  used  for  regulariza¬ 
tion.  CTRs  of  quadratic  signals  calculated  from  various  or¬ 
ders  of  TS\  D  solutions  are. collected.  With  a  defined  range 
of  system  orders,  a  CTR  plane  as  a  function  of  truncation 
parameters  and  system  orders  is  determined.  Filter  coef¬ 
ficients  lor  the  quadratic  imaging  generation  are  obtained 
from  a  truncation  parameter  and  a  system  order  that  give 
the  highest  CTR  value  in  the  CTR  plane.  Of  course,  the 
rs\  D  approach  can  be  used  to  obtain  t  he  coefficients  of 
the  quadratic  kernel.  HqU, ji),  -for  a  predetermined  filter 
order  ?n.  This  may  be  necessary  from  the  implementation 
point  of  view,  when  the  size  of  the  kernel  is  to  be  kept  at 
manageable  level.  The  results  presented  in  this  paper  are 
obtained  with  low-order  filter  to  emphasize  the  practicality 
of  the  Volterra  niter  approach. 

The  quadratic  image  is  produced  by  applying  the 
quadratic  fiber  coefficients  to  the  beamformed- RF- data 
throughout  the  standard  B-rnode  image  to  estimate  the 


f\nd 

Fig  2  A  flowchart  of  the  aigoriihm  for  quadratic  image  generation, 
quadratic  component 

m-lm  — 1 

uq(ti  4  1)  =  T.  V  «("  -j)u{v  -  I,  )Lq{].  /.•). 

3=0  k=j  ! 15 J 

where  liQij,  k)  is  the  estimated  quadratic  kernel  (extracted 
from  Ii.mnls)-  A  flowchart  of  this  algorithm  is  shown  in 

Fig-  2.  . 

V- 

D.  Harmonic  Imaging 

Second  harmonic  imaging  is  based  on  filtering  t-lie  RF 
data  with  a  zero-phase  linear  BPF  centered  ar  twice  the 
fundamental  with  restricted  bandwidth  to  minimize  the 
overlap  between  the  fundamental  and  the  second  liai- 
moiiic.  The  approach  is  appropriate  for  native  harmonic*. 


rm.KPATTARANONT  AND  EBBIM:  POSTBEAMFORMINO  SECOND-ORDER 


VOLTERRA  FILTER 


imaging.  However.  When  imaging  UCA.  the  echo  data 
.rom  the  UCA  regions  tends  to  have  wider  bandwidth 
when  compared  with  echoes  from  tissue  regions.  There¬ 
fore.  to  obtain  the  highest  possible  contrast  with  a  lin¬ 
ear  harmonic  filter,  we  varied  both  the  center  frequency 
and  the  bandwidth  of  the  filter  and  computed  a  CTR 
plane  as  a  function  of  these  parameters.  This  approach 
produced  superior  results  compared  to  standard  SH  imag¬ 
ing.  Harmonic  imaging  results  shown  in  Section  IV  are 
obtained  with  the  optimal  linear  BPF  designed  using  the 
Parks-McClellan  algorithm  [18].  All  linear  filters  used  in 
this  paper  were  implemented  in  a  zero-phase  realization 
«lb(n)  =  k(v)  *  h(-n)  *  j(n)  or  Y(f)  =  (/)|2A'(/))  The 

order  of  the  filter  was  chosen  to  achieve  50  dB  stopband 
attenuation  and  O.o  passband  ripple. 
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Fig.  3.  Schematic  for  the  imaging  .setup  for  the 
surrounded  by  the  UCA. 
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HI  Materials  and  Methods 

Imaging  results  shown  below  demonstrate  the  contrast 
enhancement  achieved  by  the  post beamfor m in g  SVF  in  the 
context  of  contrast- agent  imaging  in  tissue-mimicking  me¬ 
dia.  Two  imaging  targets  with  two  imaging  probes  are  used 
to  demonstrate  generality  and  robustness  of  the  approach. 

,-t  Contrast  Agent 

^  The  contrast  agent.  BR14  (Bracco  Research  S.A.. 
Geneva.  Switzerland),  was  used.  The  BR14  is  a  new  ex¬ 
perimental  agent  that  consists  of  high  molecular  weight 
gas  bubbles  encapsulated  by  a  flexible  phospholipid  shell 
A  0.125  mL  sample  of  BR14,  prepared  with  5  mL  of  0.9% 
saline,  was  diluted  in  500  mL  of  0.9%  saline  leading  to 
a  1:4000  dilution.  In  addition  to  BR14.  cellulose  particles 
(Sigma  Cell  Type  20,  Sigma  Chemical  Co..  St.  Louis,  MO) 

^  ere  used  as  linear  scatterers  in  the  flow  channel  target  de¬ 
scribed  below. 

B.  Imaging  Targets 


.Fig.  4  The  imaging  set.up  for  the  flow  phantom. 


small  backscatter  component  present,  due  to  the  beam  re¬ 
flections  at  the  beaker,  but  this  is  hard  to  quantify  in  the 
presence  of  the  contrast  agent..  Despite  this  limitation,  it 
is  interesting  to  compare  the  nature  of  the  image  pixels  in 
this  region  from  the  three  imaging  methods. 


Images  of  two  targets  from  two  experimental  setups 
were  used  to  evaluate  the  performance  of  the  quadratic 
imaging  technique.  The  first  experimental  setup  is  shown 
in  Fig.  3.  The  target  is  the  L-shaped  tissue-mimicking 
phantom  in  a  beaker  containing  a  dilution  of  BR14. 
The  BR14  was  constantly  stirred  bv  a  magnetic  stir¬ 
rer.  A  MEGAS  scanner  (ESAOTE  S.p.A..  Genova,  Italy) 
equipped  with  a  2-MHz  phased  array  probe  was  used  for 
image  acquisition.  The  RF  data  acquisition  was  performed 
with  15-bit  resolution  at  40-MHz  sampling  frequency  and 
without  time  gain  control  (TGC).  This  experiment  allowed 
us  to  compare  echoes  from  three  different  locations:  the 
L-shaped  tissue-mimicking  phantom,  the  BR14.  and  the 
echo- free  region  visible  in  the  scan.  The  significance  of  the 
latter  is  that  any  signal  components  observed  in  this  region 
are  largely  artifact?  from  beamfonning  and/or  reverbera¬ 
tion.  Jr  should  be  noted,  however,  that  there  may  be  a 


The  second  experimental  setup,  shown  in  Fig.  4  was 
used  in  obtaining  images  of  a  flow  phantom  (Model  524: 
ATS  Laboratories,  Inc..  Bridgeport.  CTj  containing  four 
flow  channels  with  diameters  2,  4.  G.  and  S  mm  embedded 
in  rubber-based  tissue  mimicking  material.  The  flow  phan¬ 
tom  was  connected  to  a  flow  system  with  a  roller  pump 
(Model  77200-G0:  Cole-Paruier  Instrument  Co  .  Vernon 
Hills.  IL).  Subsequently,  the  diluted  8R14  and  cellulose 
particles  were  circulated.  In  addition,  the  diluted  BR14 
was  constantly  stirred  in  a  beaker  using  a  magnetic  hot 
plate  stirrer  (EW-84303-*20:  Corning  Inc  Corning.  NY). 

This  experiment  was  designed  to  compare  linear 
backscottered  signals  from  cellulose  in  the  8- trim  diame¬ 
ter  flow  channel  with  nonlinear  backseat  ter  ed  signals  from 
BR14  in  the  6-nun  diameter  flow  channel.  The  RF  data 
were  recorded  and  saved  for  later  processing  by  the  Tech- 
nos ‘MP-  ultrasound  system  s ESAOTE  S.p.A..  Genova. 
Italy)  with  a  convex  array  prob-»  iOA42L  ESAOTE  S.p.A.. 
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Genova.  Italy)  located '  perpendicularly  to  both  the  (JCA 
and  the  cellulose  flow  channels.  The  RF  data  were  acquired 
with  16-bit  resolution  at  20-MHz  sampling  frequency  and 
without  rCiC  A  3- cycle  pulse  at  3.13  MHz  was  transmit¬ 
ted  to  produce  standard  B-mode  images. 


C.  Quantitative  Analysis 

After  the  RF  data  corresponding  to  standard  B-mode 
images  were  collected  from  experimental  setups  described 
above,  quadratic;  images  were  obtained  using  the  algorithm 
described  in  Section  II-C.  For  comparison,  harmonic  im¬ 
ages  were  produced  by  filtering  the  standard  B-mode  im¬ 
age  with  bandpass  filters  with  center  frequencies  and  frac¬ 
tional  band  widths  chosen  to  maximize  the  CTR  (defined 
by  (11))  as  described  in  Section  1I-D. 

1.  CTR  Computations:  The  CTR  was  used  for  com¬ 
parison  between  the  image  data  in  the  RF  domain  before 
scan  conversion  as  described  in  (11).  Size  and  location  of 
reference  regions  for  each  target,  are  as  follows.  For  the 
L-shape  phantom.  45-mm  axial  segments  of  20  adjacent 
A-lines  on  the  tissue  side  and  the  same  number  on  the 
contrast  side  were  used.  This  phantom  provided  an  excel¬ 
lent  opportunity  to  compute  CTR  at  multiple  regions  in 
rhe  contrast  and  tissue  regions  with  the  same  depth  and 
rhe  same  beam  angle  from  the  axis  of  the  imaging  array. 
The  calculation  region  extends  by  45  mm  in  the  axial  di¬ 
rection  (  from  50  mm  to  9o  mm)  and  comprises  20  adjacent- 
A -lines  of  RF  data  before  scan  conversion.  The  CTR  val¬ 
ues  <u e  calculated  from  cells  with  3.8-mm  axial  extent  with 
30  connected  A- lines  and  b0%  overlap  in  hath  directions. 
For  tlie  flow  phantom,  one  contrast  region  is  in  UCA  flow 
channel  and  the  tissue  region  located  between  two  flow 
channels  with  the  same  axial  extent  (3  nun  in  this  case). 
Twenty  A-line  segments  from  the  each  region  were* used. 


-4mi„  and  .4max  are  minimum  and  maximum  values,  re¬ 
spectively.  Both  .4nm,  and  .4niax  are  carefully  computed 
so  that  histograms  of  images  cover  256  gray  levels  without 
saturation.  This  was  achieved  bv  finding  the  average  of  the 
minima  (maxima)  of  all  image  lines  after  median  filtering 
This  was  done  t.o  ensure  that  the  display  range  is  not  s<%\ 
by  extreme  values  of  .4mir,  or 

From  gray-level  images,  the  contrast  uuio  iCH)  :15! 
used  as  the  contrast  measurement  between  any  two  regions 
is  given  by 


/‘2  .  *>  *  *  *  * 

V  &  ]  “t" 

where  I\  and  U  are  the  average  of  gray  levels,  and  and 
<72  are  the  corresponding  standard  deviations  in  the  first 
region  and  the  second  region,  respectively 

3.  Histogram  and  Receiver  Operating  ChuracU  rist/c 
Analysis:  In  the  imaging  results  shown  below,  we  evaluate 
the  CTR  and  CR  based  on  regions  in  the  image  represen¬ 
tative  of  the  UCA  and  tissue  regions  with  the  same  num¬ 
ber  of  pixels  and  at  the  same  depth.  When  the  number  of 
pixels  in  these  regions  is  sufficiently  high  t-o  produce  mean¬ 
ingful  statistics,  histograms  of  the  pixel  data  aie  produced 
to  demonstrate  these  statistics.  In  addition,  mvivei  oper¬ 
ating  characteristics  (ROC)  analysis  [16]  is  performed  as  » 
simple  method  for  classificat  ion  of  the  different  regions.  In 
addition,  when  additional  regions  can  be  identified  teg. 
low-scattering  regions)  in  any  image,  the  CR  between  tis¬ 
sue  and  such  regions  is  also  evaluated.  Please  note  that  the 
ROC  curves  are  included  in  this  paper  to  further  quantify 
the  degree  of  overlap  between  different  histograms  from 
tissue,  contrast,  and  low-scattering  regions. 


d  ('<>v frast  Ratio:  As  can  be  seen  from  results  given 
below,  the  dynamic  range  (in  decibels)  of  images  from 
quadratic  components  is  approximately  twice  the  dynamic 
range1  of  standard  B-morie  and  SH  images.  In  order  to 
account  for  image  perception  on  standard  $-lut  display, 
all  images  are  represented  with  their  full  dynamic  range 
mapped  to  256  gray  levels.  That  is.. 


where  A  is  the  log  magnitude  pixel  values,  given  by 


D.  Resolution  Measurements 


We  used  correlation  lengths  calculated  from  the  2D  au¬ 
tocorrelation  of  echoes  in  tissue  regions  as  described  in  jl7l 
to  measure  spatial  resolution.  Intensity  images  are  scan 
converted  and  uniform  speckle  regions  were  identified  to 
compute  the  average  speckle  correlation  cell  size 


Cj{x\y) 

0(0.0) 


eh-  dy. 


(10) 


where  Tilx)  denotes  the  Hilbert  transform  of  data  vector, 
r.  and  j-TCr{  >  0)  is  some  constant  (typically  maximum  am¬ 
plitude).  The  data  vector  is  the  digitized  RF  for  standard 
B-mode  and  rhe  oiuput  of  the  harmonic  or  the  quadratic 
filter  for  the  harmonic  and  quadratic*  images,  respectively. 


where  Cj(x.y)  is  the  2B  correlation  function  of  the  inten¬ 
sity  autocovariance  function,  and  A*  and  )'  are  taken  ro 
be  sufficiently  large  t,o  allow  the.  magnitude  of  the  auto 
covariance  to  drop  to  negligible  levels.  The  vertical  and 
horizontal  correlation  cell  sizes,  representing  the  axial  and 
lateral  resolution,  respectively,  were  used  to  compare  spa¬ 
tial  resolution  for  the  throe  imaging  techniques 
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IV.  Results 

.*1.  L  -Shaped  Phantom 

Fi».  a  shows  tlie  average  spectra  of  typical  echoes  from 
the  contrast  and  tissue  regions  of  the  L-shaped  phantom. 
The  average  spectra  were  calculated  by  averaging  win¬ 
dowed  periodogram  of  every  echo  line  in  regions  described 
in  Section  HI-C.l.  These  spectra  suggest  that  eclioes  from 
tlie  contrast  region  have  broader  bandwidth  compared  to 
the  echoes  from  tissuelike  region.  As  mentioned  in  Sec¬ 
tion  11-D.  this  suggests  that  the  best  linear  filter  for  con- 
tnist.  enhancement  may  not  be  the  standard  SIT  filter. 
Rather,  a  genera]  BPF  with  center  frequency  and  relative 
bandw  idth  designed  to  maximize  the  CTR  may  be  sought. 
Based  on  this  approach,  we  have  designed  a  zero-phase 
BPF  with  center  frequency  of  3.2  MHz  and  relative  band¬ 
width  of  30%  to  demonstrate  the  performance  of  linear 
post  beamforming  harmonic  imaging.  The  quadratic  filter 
was  derived  from  echoes  from  the  tissue  region  according  to 
the  algorithm  described  in  Section  II-C  (with  P  =  Q  =TiO 
in  Fig.  2}.  Panels  (a),  (b).  and  (c)  of  Fig.  G  show  images  ob¬ 
tained  using  standard  B-inode,  linear  BPF.  and  quadratic 
imaging  methods,  respectively.  Due  to  differences  in  dy¬ 
namic  ranges  for  the  three  methods,  each  image  is  dis¬ 
played  with  its  full  dynamic  range 'as  can  be  seen  from 
the  decibel-level  scale  bars  shown  in  accordance  with  (1G). 
fig.  (5(d)  shows  fou^  rectangular  regions  used  for  charac¬ 
terization  of  the  imaging  results.  Regions  „4i  and  A2  are 
representative  of  tissue  and  contrast  echoes,  respectively. 
Regions  Bi  and  B2  *are  representative  of  tissue  and  echo- 
free  region,  respectively.  These  regions  are  used  for  CTR 
calculations  as  well  as  evaluating  histograms  of  the  recon¬ 
structed  images.  Due  to  the  structure  of  the  target,  regions 
B:  and  B2  could  not.  bo  placed  at  the  same  depth  in  the  im¬ 
age.  Nevertheless,  region  Bi  is  representative  of  tissue  re¬ 


sponse  and  provides  a  valid  baseline  for  comparison  of  con 
trast  with  the  echo-free  region.  Bo.  The  standard  B-mode 
image  (Fig.  6(a))  was  obtained  from  the  digitized  RF  echo 
data  and  displayed  with  a  dynamic  range  of  60  dB  The 
structure  of  the  L-shaped  phantom  can  be  recognized  with 
low  contrast  between  the  UCA  and  the  1  issue-mimicking 
regions.  The  strong  specular  reflection  at  the  top  of  rhe 
UCA  region  is  due  to  a  boundary  Inver  formed  In  the 
agent  at  the  interface  with  the  tissue-mimicking  medium 
A  careful  examination  of  the  image  reveals  a  higher  level 
of  contrast  between  the  two  regions  at  dose  range.  A  CTR 
value  of  6.69  dB  was  computed  based  on  echo  signals  from 
regions  A2  and  .4]  (Fig.  6(d)). 

The  linear  harmonic  BPF  image  shown  m  t  ig.  (i(b)  was 
produced  by  filtering  the  RF  data  with  bandpass  filter  cen¬ 
tered  around  the  second  harmonic  (3  2  MHz)  with  a  frac¬ 
tional  bandwidth  of  30%.  The  center  frequency  and  the 
fractional  bandwidth  of  the  harmonic  filter  were  chosen  to 
optimize  the  CTR  value  for  regions  .4]  and  A>  in  Fig  Old: 
The  CTR  value  for  the  harmonic  image  (18.2  dB)  is  con¬ 
sistent.  with  the  perceived  enhancement  clearly  visible  in 
the  image.  One  can  also  observe  the  loss  of  the  specular 
reflections  at  the  top  right  and  bottom  left  of  the  harmonic 
image  This  is  typical  since  these  echoes  have  significant 
low  frequency  components.  As  expected,  the  speckle  m  the 
tissue  region  appears  finer  than  that  of  the  standard  B 
mode. 

The  quadratic  image  obtained  using  the  algorithm  de¬ 
scribed  in  Section  II-C  is  shown  in  Fig.  6(c)  The  CTR 
for  this  image  is  21.3  dB  indicating  contrast  enhancement 
over  both  standard  B-inode  and  harmonic  images.  One 
important  feature  of  the  quadratic  image  is  t  he  increased 
dynamic  range  compared  with  3-mode  and  harmonic  im¬ 
ages  This  increased  dynamic  range  results  in  contrast  en¬ 
hancement  without  loss  in  image  features,  such  as  specular 
reflections,  which  may  be  of  diagnostic  value  m  some  cases. 

Another  interesting  comparison  among  the  three  ini 
age  results  shown  in  Fig.  G(a)-(c)  is  the  echo-free  region 
visible  at  the  right  edge  of  the  scan  One  can  see  that 
contrast  between  tissue  and  this  region  in  the  quadratic 
image  is  superior  to  that  from  the  B-mode  image,  vheieas 
this  contrast  in  the  harmonic  image  is  the  lowest.  This 
result  is  significant  because  echo  signals  from  tins  region 
are  largely  artifacts  (due  to  beamformmg  and  /ot  rewi  l  or¬ 
ations).  Quantitative  measures  of  this  cont  rast  enhance¬ 
ment  betw-een  the  tissue  and  echo- free  regions  are  given 
below. 

To  further  illustrate  the  imaging  results  given  in  Fig.  6. 
vertical  and  horizontal  lines  through  images  are  plotted 
in  Fig.  7.  Axial  lines  from  the  three  imaging  techniques 
through  the  UCA  region  are  shown  in  Fig.  7(a).  One  can 
see  that  the  specular  reflector  is  all  but  eliminated-  in  riic 
harmonic  image,  while  it  remains  visible  for  the  other  two 
methods.  Further,  the  axial  line  from  the  quadratic  image- 
shows  no  apparent,  loss  in  the  axial  resolution  at  the  spec¬ 
ular  reflection  location.  Lateral  lines  at  rhe  90- mm  depth 
from  the  three  imaging  techniques  are  shown  in  Fig  7  b  . 
These  lateral  lines  pass  through  the  UCA.  tissue,  and  :he  * 
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Ki?.  0.  im»8«  of  the  beaker  containing  the  tissue-mimicking  L-sliaped  phantom  surrounded  bv  the  UCA  from  three  imamn- 
(a)  Standard  B-mode.  (b)  Linear  harmonic  BPF.  (c)  Quadratic,  (d)  Box<*  indicate  regions  for  CTR  and  histogram  calculations 
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echo-free  regions.  One  can  see  the  contrast  enhancement 
from  t  he  lateral  line  of  the  harmonic  image  compared  with 
1  he  lateral  line  of  the  standard  B-mode  image.  Most  signif¬ 
icantly.  however,  one  can  see  that  the  quadratic  image  pro¬ 
vides  the  maximum  separation  between  the  three  regions 
without  apparent  loss  of  lateral  resolution.  This  was  con¬ 
firmed  by  computing  the  correlation  cell  size  in  the  scan- 
converted  intensity  images  according  to  (19)  [17]  from  the 
Three  imaging  methods  in  the  tissue  region  (59  to  73  mm 
axial  and  —20  to  —14  mm  lateral  in  Fig.  6).  These  values 
are  reported  in  Table  I. 

B.  Flow  Channel  Phantom 

Hg.  S  shows  the  average  spectra  of  typical  echoes  from 
tlm  contras?  and  tissue  regions  of  the  flow  phantom  de¬ 
scribed  m  Section  IIJ-B.  The  average  spectra  are  calcu¬ 
lated  by  averaging  windowed  periodogram  of  even*  echo 


line  in  regions  described  in  Section  III-C.I.  As  with  the 
L-sliaped  phantom  result,  the  echoes  from  the  UCA  re¬ 
gion  exhibit  broader  bandwidth  than  those  from  tissue  re¬ 
gion.  The  standard  B-mode  image  of  the  flow  phantom 
consisting  of  the  UCA  and  cellulose  in  flow  channels  was 
acquired  using  the  experimental  setup  described  in  Sec 
rion  I1I-B.  The  image  is  shown  in  Fig.  9(a).  One  can  sot* 
the  baekscattered  enhancement  due  to  the  UCA  in  the  (>• 
mm  flow  channel  (55-  to  Gl-mm  axial  and  -6  to  0  nun 
lateral)  compared  with  the  surrounding  tissue  with  CTR 
9.8  dB.  The  CTR  was^calculated  based  on  echo  signals 
from  the  6-mm  flow*  channel  excluding  the  specular  reflec¬ 
tors  (Box  Ci)  and  echoes  from  tissue  at  the  same  depth 
(Box  Co).  On  the  other  hand,  the  echoes  from  the  cellulose 
in  the  8- mm  flow  channel  (53  to  61  rrim  axial  and  17  n. 
25  mm  lateral)  are  weaker  than  those  from  the  surrounding 
tissue.  A  linear  scar  terer- to- tissue  ratio  (LTR)  was  calcu¬ 
lated  nased  on  echo  signals  from  the  S-mni  now  chaum-l 
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TABLE! 

C  h.  Cl  R.  AND  Correlation-  Cell  Size  Vam  es  for  the  L- Shared  Imagine  Target  The  CTR  Valves  auk  Given  as  Mean  r. 

Decibels  with  Standard  Deviation  in  Parentheses. 
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I-  ig  <  Lines  through  images  in  Fig.  G  from  t hree  imaging  techniques 
Thin-  B-mode.  Dash:  Harmonic.  Thick:  Quadratic.  (Top)  Axial  lines 
through  the  UOA  region  (Bottom)  Lateral  lines  (at  depth  90  mm) 
through  the  tissue-mimicking.  UCA.  and  echo-free  regions. 


(Box  Cs)  and  echoes  from  tissue  at  the  same  depth  iBo\ 
C2).  The  LTR  was  -8.8  dB,  calc  ulated  based  on  '  II »  fm 
8-iuui  channel. 

The  harmonic  image  shown  in  Fig  9(b)  was  produced 
by  filtering  the  RF  data  with  bandpass  filter  centered 
around  4.8  MHz  with  a  fractional  bandwidth  of  30 %  The 
center  frequency  and  the  fractional  bandwidth  of  the  hat 
monic  niter  were  chosen  to  optimize  the  CTR  value  usimi 
the  method  described  in  Section  II -D  for  the  UCA  and 
tissue  regions  shown  in  Fig.  9(d).  The  CTR  value  for  rht 
harmonic  image  is  15  dB.  which  is  consistent  with  the  per 
oeived  enhancement  clearly  visible  in  the  displayed  image, 
it  is  interesting  to  note  the  loss  of  the  specular  reflections 
from  the  contrast  flow  channel  but  not  from  .the  lineai  scat 
terer  channel.  It.  is  also  interesting  to  note  that  the  speckle 
in  the  tissue  region  is  finer  than  that  of  the  standard  B- 
.  niodc.  This  indicates  that  the  optimization  of  the  harmonic 
filter  resulted  in  a  larger  fractional  bandwidth  than  that 
of  the  fundamental  component.  The  LTR  for  the  harmonic- 
component  is  now  -7.6  dB.  That  is.  the  harmonic  image 
does  not  preserve  the  low-scattering  region  as  it  improves 
the  contrast  between  the  UCA  and  tissue  regions 

Fig.  9(c)  shows  the  quadratic  image  obtained  using  the 
algorithm  described  in  Section  II-C  (with  P  =  Q  =  12 
in  Fig.  2).  Compared  with  the  standard  B-mode.  the  en¬ 
hanced  visualization  of  the  UCA  in  the  small  flow  channel 
surrounded  by  the  tissue  can  be  seen.  This  enhanced  visu¬ 
alization  is  also  consistent  with  the  higher  CTR  (22.2  dB). 
In  addition,  one  can  see  the  reduction  in  the  echogenic¬ 
ity  of  the  cellulose  in  the  8- mm  flow  channel.  This  is  also 
reflected  in  the  value  of  LTR  (-15.4  dB)  computed  foi 
this  channel.  This  result  demonstrates  the  ability  of  the 
quadratic  filter  to  separate  the  tissue  from  both  the  UCA 
and  the  low- scattering  cellulose  region  simultaneously.  We 
can  also  see  the  preservation  of  specular  reflections  from 
both  the  UCA  and  cellulose  flow  channels.  This  can  be  seen 
more  quantitatively  in  Fig.  10.  which  show*?  the  decibel 
values  for  axial  and  lateral  lines  through  images  in  Fig  -9. 
Axial  lines  through  the  center  of  the  cellulose  and  the  UCA 
flow-  channels  are  shown  in  Fig.  10(a)  and  Fig  10(b).  re¬ 
spectively.  Lateral  lines  Through  the  center  of  both  the  cel¬ 
lulose  and  the  UCA  flow  channels  are  shown  in  Fig.  10(c). 
One  can  see  the  contrast  enhancement  m  different  *  me¬ 
dia  (he.,  the  UCA.  the  cellulose,  and  the  tissue i  from  the 
quadratic  image  without-  loss  in  spatial  resolution  In  par¬ 
ticular.  for  the  B-mode  image,  the  signal  from  the  con¬ 
trast  channel  is  13  dB  above  the  tissue  ievei  and  3! 1  -iB  ■ 
above  the  cellulose  channel.  On  the  other  hand,  foi  \u*> 
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B-mode 
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Quadratic 
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The  defined  regions 
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Hg  *  Images  cf  the  flow  phantom  using  (a)  Standard  B-mode.  (b)  Harmonic,  (c)  Quadratic,  (d)  Boxes  C,  .To.  and  O*  indicate  regions  for 
CTK  *nd  CR  calculat,ons*  Note  U,e  hvpcrechoic  region  from  the  UCA  in  the  6-mm  flow  channel  (55  mm  to  61  mm  axial  and  -6  trim  to  0 
lateral :  and  the  bypoechotc  region  from  the  cellulose  (53  mm  to  01  mm  axial  and  IT  mm  to  25  mm  lateral). 


quadratic  image,  the  signal  from  the  contrast  channel  is 
30  clB  higher  than  tissue  and  67  aB  higher  than  the  cel¬ 
lulose  channel.  Finally,  examination  of  the  strong  specular 
reflections  along  the  axial  lines  shown  in  Fig.  10  suggests 
that  the  quadratic  filter  preserves  the  resolution  of  the 
system.  This  was  confirmed  by  computing  the  correlation 
cell  size  in  the  scan-converted  intensity  images  according 
to  (19)  [17]  from  the  three  imaging  methods  in  the  tissue 
region  (52  to  65  nun  axial  and  4  to  17  lateral  in  Fig.  7), 
These  values  are  reported  in  Table  II. 

C.  Classification  Results' 

To  give  the  reader  a  quantitative  idea  of  the  ability  of 
(lie  three  different  imaging  methods  to  separate  the  dif¬ 
ferent  regions  in  the  imaging  targets,  we  use  histograms 
from  representative  areas  of  these  regions.  For  example. 


for  the  tissue-mimicking  L-shaped  phantom  surrounded 
by  the  UCA  in  the  beaker,  regions  from  different  media 
are  defined  in  Fig.  6(d)  as  follows:  the  region  A1  and  A2 
(the  tissue  and  UCA  regions,  respectively)  and  the  re¬ 
gion  Bl  and  B2  (the  tissue  and  echo-free  regions*  respec¬ 
tively).  After  images  from  three  different  imaging  tech¬ 
niques  are  scan  converted  and  represented  with  256  levels 
of  gray  covering  their  full  dynamic  range,  the  histogram 
of  each  region  is  determined  and  shown  in  Fig  11.  His¬ 
tograms  from  regions  A-l  and  A2  are  shown  in  Fig.  1] 
(left-hand  side).  One  can  see  the  degree  of  overlap  be- ' 
tween  the  histograms  is  highest  for  the  standard  B-mode 
image.  On  the  other  liand.  the  corresponding  histograms 
are  well  separated  for  the  harmonic  and  quadratic  images. 
This  is  consistent  with  the  increased  level  of  contrast  per¬ 
ceived  .from  direct  visualization  of  the  images  in  Fig.  6. 
Histograms  produced  from  regions  ‘**1  and  B2  an*  shown  in 
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Fig.  10.  Lines  ihrougli  the  standard  B-mode  (thin),  harmonic  (dashed),  and  the  quadratic  (thick)  from  images  shown  in  Fig.  9  ;Toj>;  Axial' 
lines  through  the  center  of  the  cellulose  flow  channel.  (Middle)  Axial  lines  through  the  center  of  the  UCA  flow  channel.  (Bottom)  Lateral 
iin-s  through  the  center  of  both  two  flow  channels 


TABLE  I] 

(.  H.  Cm.  AM)  C ’OK RELATION  C’ELL  SIZE  VALUES  FOR  THE  Fl.OW  CHANNEL  TARGET.  THE  CTR  VALVES  ARE  GIVEN  IN  Dt/  IREl. 


Imaging  method 

UCA/Tissue 

Tissue/ Cellulose 

CTR 

Sz 

ST 

15-mode 

3 .57 

1.33 

9.8 

1.0 

1.2 

SH 

2.12 

1.21 

35.0 

0.6 

0.9 

Quadratic 

2  46 

1.65 

22.2 

0.9 

1.1 

Pig.  J 1  (right-hand  side).  Grav-level  histograms  of  regions 
Bl  and  B2  produced  from  harmonic  image  has  higher  de¬ 
gree  of  overlap  than  those  from  the  standard  B-mode  and 
the  quadratic  images.  This  is  further  quantified  by  ROC 
cuives  [3Gj  obtained  from  histograms  of  region  Al  (Tis¬ 
sue)  and  A2  (UCA)  and  region  Bl  (tissue;  and  B2  (echo- 
free)  are  shown  in  Fig.  12(a)  and  Fig.  12(b).  respectively. 
As  can  be  seen  in  Fig.  12(a).  the  .4-  values  between  re¬ 
gion  Al  and  A2  from  the  standard  B-mode.  harmonic  and 
quadratic  images  are  0.8269.  0.9936.  and  0.9876.  respec¬ 
tively.  These  values  demonstrate  improved  classification 
performance  between  region  Al  and  region  A2  from  har¬ 
monic  and  .quadratic  images  over  the  standard  B-mode  im¬ 
age.  However,  as  shown  in  Fig.  12U>).  the  .4.  value  of  the 
quadratic  image  (0.9535)  presents  the  best  classification 
performance  of  tissue  and  echo-free  regions  among  three 
imaging  Techniques,  while  tin*,  classification  performance 
from  the  harmonic  image  (A:  =  0.8724}  is  inferior  ro  that 


from  the  standard  B-mode  image  (Az  ~  0.9302).  These 
measurements  show  that,  the  quadratic  image  provides  int 
proved  separation  of  the  UCA  from  the  tissue  mimic  com¬ 
parable  with  SH  performance.  At  the  same  time,  signifi¬ 
cant  improvement  in  separation  between  tissue  mimic  and 
echo- free  regions  is  achievable  in  the  quadratic  image  over 
both  the  standard  B-mode  and  SH  images. 

Contrast  ratios  determined  from  corresponding  gray 
scale  images  in  Fig.  6  are  shown  in  Table  3.  For  contrast- 
comparison  between  the  UCA  and  tissue  regions.  CRs  ob¬ 
tained  from  regions  A! "and  A2  demonstrate-  the  contrast 
enhancement  of  the  harmonic  (2.56)  and  the  quadratic 
(2.16)  images  over  the  standard  B-mode  image  0  91).  On 
ike  other  hand,  one  can  see  that  the  CR  between  tissue  and 
echo-free  regions  of  the  harmonic  image  si  12)  is  iuferu*?* 
to  that  of  the  standard  B-mode  image  (1  51 ).  However,  the 
quadratic  image  gives  the  maximal  CR  ‘(1.66i.  These  CR 
values  agree  with  both  the  visualization  of  images  shown 
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Bmode 


f tg.  1 1  Gray -level  histograms  produced  from  images  shewn  in  Fig. 
“ti  the  left-hand  side  are  produced  from  region  A I  ?uid  A2  Thick: 
right -hand  side  are  produced  from  region  Bl  and  B2.  Thick:  Region  ! 


Bmode 


Gray  level 

6.  Top:  B-rnode.  Middle:  Harmonic.  Bottom:  Quadratic.  Histogiam- 
legion  A1  (Tissue)  Thin:  Region  A2  (Contrast)  Histograms  on  tie- 
II  (Tissue).  Thin:  Region  B2  (Air). 


I  ig  12  Ihe  ROC  curves  produced  from  gray-level  histograms  in  Fig.  11  of  three  imaging  techniques.  Thick:  B-inodt*.  Da*»h:  JlarintMiir 
1 !;:::  Quadratic.  The  ROC  curve  obtained  from  the  UCA  find  tissue  regions  is  shown  on  the  left-hand  side.  The  ROC  t  urve  obtained  from 
•  he  vctL-rei  and  tissue  regions  is  shown  on  the  tight  -hand  side.  .  • 
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in  Pig.  6  and  the  degree  of  overlap  between  the  histograms 
shown  in  Fig.  11.  These  results  have  demonstrated  that, 
the  quadratic  image  provides  the  perceived  contrast  en¬ 
hancement  not  only  between  the  UCA  and  tissue  regions 
but  also  between  the  tissue  and  the  echo-free  regions.  We 
note  here  that  similar  analysis  was  performed  for  the  flow 
channel  data  and  that  the  results  are  in  full  agreement 
with  the  results  shown  in  this  section.  The  CR.  CTR.  and 
correlation  cell  si2e  values  for  both  the  L-shaped  and  How 
phantoms  are  summarized  in  Tables  1  and  2. 

Similarly,  images  of  the  UCA  and  cellulose  How  channels 
m  the  flow  phantom  as  shown  in  Fig.  9  are  scaled  into  8-bit 
gray -level  images.  Gray- level  histograms  between  the  UCA 
and  tissue  regions  and  between  the  cellulose  and  tissue 
regions  are  produced  and  the  corresponding  ROC  curves 
arc  determined.  Although  the  number  of  data  from  each 
region  is  quite  small  for  meaningful  statistical  analysis, 
results  from  both  CRs  and  Az  show  tendency  that  the 
quadratic  image  enhances  the  contrast  and  provides  better 
classification  performance  compared  with  the  standard  B- 
modo  image.  Table  II  gives  the  CR  values  resulting  from 
the  flow  channel  experiment. 

V.  Discussion 

We  have  introduced  a  nonlinear  post  beamforming,  fil¬ 
tering  algorithm  for  ultrasonic  pulse-echo  imaging  based 
on  the  Volte?  ra  filter  model.  The  main  goal  of  this  paper 
to  introduce  the  mathematical  basis  for  deriving  the 
filter  coefficients  from  standard  beam  formed  RF  data  ob¬ 
tained  by  commercial  scanners.  In  addition,  we  have  pre¬ 
sented  imaging  results  from  two  laboratory  contrast  tar¬ 
gets  to  illustrate  the  nature  of  the  gray  scale  images  ob¬ 
tained  with  the  quadratic  component  compared  with  stair 
‘laid  echo  signals  and  harmonic  images.  Images  from  the 
quadratic  components  were  obtained  by  applying  a  sin¬ 
gle  filter  derived  from  echoes  from  the  contrast  region.  We 
have  decided  to  do  this  to  demonstrate  that  the  method  is 
quite  lobust  in  the  sense  that  the  derived  filter  is  applied 
Throughout  the  image  to  produce  images  free  of  artificial 
inliomogeneity.  This  is  a  desirable  method  from  the  im- 
plrmentation  point  of  view,  especially  for  real-time  imple¬ 
mentation. 

For  a  fair  c  omparison,  all  images  presented  in  this  paper 
were  normalized  to  their  full  dynamic  range  and  displayed 
using  256  levels  of  gray.  Classification  results  based  on 
histogram  characterization  of  contrast  and  tissue  regions 
show  that  the  quadratic  images  produce  nearly  twofold 
increase  in  CR  values  (compared  with  standard  echo  im¬ 
ages.)  without  loss  of  image  features  (compared  to  har¬ 
monic  images).  Classifying  the  tissue  and  echo-free  regions 
demonstrated  the  improved  performance  with  respect  to 
harmonic  imaging.  This  is  probably  due  the  vulnerability 
of  harmonic  images  to ’noise  and  beainformmg  artifact  s, 
especially  when  the  sidclobe*  of  the  transmit  beam  are  in 
contrast  regions. 

though  we  have  focused  on  the  comparison  be-' 
rwi-cn  gray-scale  images,  there  are  some  interesting  prop¬ 


erties  of  the  quadratic  signal  components  that  mar  reveal 
important  information  on  the  nature  of  the  objects  pro¬ 
ducing  the  echo  signals.  For  example,  the  frequency  com¬ 
ponent  at  /  in  the  quadratic  signal  component  is  a  result 
of  all  frequency  components  f\  and  fn  from  tin*  echo  Mg- 
nal  such  that  f\  +  h  -  /•  weighted  by  HQ(f  i.  f This 
frequency  coupling  results  in  quadratic  signal  components 
with  high  SNR  values  due  to  rejection  of  addiiiw  nc-i.-t* 
An  interesting  question  is  whefhei  the  quadratic  compo 
nents  are  more  directly  related  to  the  composition  of  ihe 
echoes  in  terms  of  coherent  and  diffuse  scattering  Tins 
may  be  quite  significant  in  improving  the  robustness  of 
motion  tracking  and  displacement  estimation  algorithms 
In  this  paper,  the  filter  is  designed  based  on  nonlin¬ 
ear  echoes  from  the  tissue-like  medium.  However,  the  til 
ter  does  not  distinguish  between  non  linearities  from  UCA 
and  those  from  tissue.  The  contrast  enhancement  is  main  I’, 
due  to  the  higher  level  of  quadratic  component  { relative  to 
the  RF  signal  level)  in  the  contrast  regions  For  example, 
in  the  L-shaped  phantom,  on  average,  the  quadratic  sig 
rial  is  30  dB  and  50  dB  below  the  RF  in  the  contrast  mid 
tissue  regions,  respectively.  It  is  interesting  to  note  that, 
for  this  target,  the  quadratic  component  from  the  mag¬ 
netic  bead  at  the  bottom  of  the  image  is  only  3  dB  be¬ 
low  the  RF  signal.  This  suggests  that  the  quadratic  filtei 
does  not  completely  reject  quadratic  echoes  from  specu¬ 
lar  reflectors,  even  though  it  may  still  discriminate  again*! 
them.  Therefore,  for  applications  in  which  the  quadrat  it 
components  from  the  contrast  agents  are  weak  or  compa¬ 
rable  to  tissue  components,  a  method  for  separating  tin- 
two  types  of  nonlinearity  is  needed.  There  are  several  pos¬ 
sible  approaches  to  the  filter  design  problem  so  that  better 
separation  between  nonlinear  echoes  from  UCA  and  tissue. 

•  Adaptive  implementation  of  the  SVF  based  on  fast 
recursive  LS  approach  [12].  This  will  allow  foi  ihe  op¬ 
timization  of  filter  coefficients  based  on  the  local  level 
of  nonlinearity. 

•  Higher  order  filters  (e.g.,  cubic)  that  may  be  more 
sensitive  to  UCA  nonlinearity  than  tissue  nonlinear¬ 
ity’.  This  is  due  to  the  observation  that,  under  nor¬ 
mal  imaging  conditions,  tissue  nonlinearity  is  at  most 
quadratic. 

•  Synthesis  of  pulse  sequences  to  excite  contrast  nii- 
crobubbles  to  maximize  their  nonlinear  response.  This 
is  motivated  by  the  recent  trend  in  conti  ast -agent, 
imaging,  which  calls  for  the  use  of  super  low  value* 
of  MI.  The  identified  quadratic  kernel  of  the  SVF  can 
be  used  for  the  design  of  these  waveforms 

We  have  chosen  to  compare  quadratic  images  with  im¬ 
ages  obtained  using  linear  harmonic  filters  rather  than  sec¬ 
ond  harmonic  filters.  This  was  based  on  the  observation 
that  the  spectra  of  echoes  from  contrast,  regions  are  typi¬ 
cally  broader  than  those  from  tissue  regions  as  can  be  seen 
from  Figs.  5  and  8.  The  application  of  a  strict  SH  filter 
based  on  the  bandwidth  from  tissue  component  typically 
produce  inferior  results  when  UCA  regions  arc  present  in 
the  imaging  field.  To  illustrate  this  point.  Fig.  13  show- 
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2nd  harmonic 


Lateral  (mm) 


Fig.  13.  LHS:  The  image  of  L-shaped  phantom  produced  using  BPF  centered  at  4  MHz  with  fractional  bandwidth  23%.  HUS:  The  image*  «M' 
how  channel  phantom  produced  using  BPF  centered  at  6  MHz  with  fractional  bandwidth  25%. 


SH  images  of  both  the  L-shaped  and  flow  phantoms  based 
on  the  spectra  from  tissue  regions.  These  images  are  .vis¬ 
ibly  inferior  to  those  obtained  using  the  linear  harmonic 
approach  described  in  Section  II-D.  We  note  that  quanti¬ 
tative  analysis  of  these  images  confirms  the  inferiority  of 
the  SH  images  shown. 


The  niters  used  for  extracting  the  quadratic-  compo¬ 
nents  for  the  two  data  sets  were  of  orders  10  and  12  for 
the  L-shaped  phantom  and  the  flow-channel  phantom,  re¬ 
spectively  In  deriving  the  filters,  we  have  assumed  that 
IiqIJJs)  =  hQ{k,j).  which  implies  that  the  filters  can  be 
implemented  with  Go  and  90  independent  coefficients,  re¬ 
spectively.  On  the  other  hand,  the  linear  bandpass  filter 
luid  77  coefficients.  This  number  is  effectively  doubled  by 
the  zero-phase  implementation  described  in  Section  II-D. 
While  it  is  not  our  objective  to  compare  the  computational 
efficiency  of  I  lie  quadrat  ic  filters,  with  the  linear  filter,  we 
note  that  the  requirements  of  the  quadratic  filters  will  not 
present  a  severe  problem  for  real-time  implementation  on 
modern  ultrasound  scanners.  We  also  note  that  efficient, 
software  and  hardware  implementations  of  the  Vol terra  fil¬ 
ter  have  been  extensively  studied  in  recent  literature  [19]. 


Finally,  even  though  the  two  imaging  targets  in  this  pa¬ 
per  contain  UCA  regions,  quadratic  images  can  be  used  in 
“native  quadratic'*  form,  in  much  the  same  way  as  SH  im¬ 
ages.  The  enhancement  of  the  contrast  between  the  tissue 
mimic  and  the  low  scattering  regions  in  both  targets  offer 
an  illustration  of  the  nature  of  quadratic  images  in  native 
mode.  We  are  currently  investigating  speckle  reduction  in 
quadratic  images  (some  speckle  reduction  can  be  observed 
in  t he  tissue  mimic  in  quadratic  images  from  both  phan-* 
tomsu  Tliis  is  the  subject  of  a  future  report. 


VI.  Conclusions 

A  nonlinear  post  beam  forming  filter  based  on  the 
Volterra  model  was  introduced  in  this  paper.  The 
quadratic  component  from  the  SVF  signal  separation 
model  was  shown  to  produce  images  with  high  contrast 
and  high  dynamic  range  without  loss  of  axial  oi  latelai 
resolution.  This  was  confirmed  by  estimating  the  correla¬ 
tion  cell  size  based  on  [17]  for  the  quadratic  images  and 
comparing  with  those  for  B-mode  and  Harmonic  images 
(Tables  I  and  II).  The  improvement  in  contrast  was  con¬ 
firmed  by  quantitative  measures,  both  on  the  RF  data  and 
the  log-compressed  image  data  after  scan  conversion.  Fur¬ 
thermore,  the  quadratic  images  were  also  shown  to  pre¬ 
serve  the  low-scattering  targets  while  improving  the  UCA 
to  tissue  contrast.  This  was  demonstrated  by  the  computed 
CR  values  (Tables  I  and  II)  and  the  Az  values  from  the 
ROC  curves  shown  in  Fig.  12. 
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